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PREFACE 
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(2)  Synthesis  and  Analysis  of  optical  waveguides  with  prescribed  TM  modes;  (3)  development 
and  testing  of  direct  scattering  solver  to  analyze  optical  waveguides;  (4)  development  of  inverse 
scattering  theory  for  the  design  of  planar  optical  waveguides  with  same  propagation  constants  for 
different  frequencies;  (5)  Analysis  of  coupling  in  multilayered  waveguides  using  inverse  scattering 
techniques;  and  (6)  Sotiton-soliton  interaction  in  nonlinear  optical  waveguides  and  bistability  in 
nonlinear  periodic  media. 

Five  journal  articles  have  been  prepared  during  the  preiod  1992-1993  based  on  the  research 
carried  out  on  this  project  and  they  are  in  different  stages  of  publication.  Two  Ph.D.  dissertations 
are  also  being  carried  out  under  this  grant.  Reprint  and  Preprint  of  selected  research  articles  are 
attached  in  the  appendix. 

The  award  of  this  grant  has  been  of  great  help  to  me  and  my  students  in  carrying  out  our 
research  and  is  gratefully  acknowledged. 
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DESCRIPTION  OF  RESEARCH  CARRIED  OUT 


1)  Synthesis  and  Analysis  of  Planar  Optical  Waveguides  with  Prescribed  TM  Modes 

An  inverse  scattering  approach  to  designing  optical  waveguides  with  prescribed  propagation 
characteristics  of  TM  modes  is  developed.  The  refractive  index  profile  of  the  waveguide  is 
formulated  as  a  solution  to  a  nonlinear  differential  equation  whose  forcing  function  is  the  potential 
obtained  from  the  application  of  the  inverse  scattering  theory.  This  method  can  reconstruct  smooth 
refractive  index  profiles  for  planar  waveguides  that  support  single  mode  or  multi  modes,  Both 
the  cases  of  zero  and  nonzero  reflection  coefficient  characterizing  the  transmission  properties  of 
waveguides  are  discussed  here.  A  direct  analysis  technique  based  on  finite  difference  scheme 
has  been  formulated  to  verify  the  results  obtained  by  inverse  scattering  method  and  they  are  in 
excellent  agreement. 

2)  Guided  Wave  Optical  Interconnects 

A  guided  wave  optical  interconnect  that  reduces  or  eliminates  clock  skew  by  ensuring  simulta¬ 
neous  delivery  of  dock  pulses  to  chips  mounted  on  a  wafer  (see  Fig.  1)  has  been  designed.  The 
interconnect  consists  of  a  multimode  planar  trunk  waveguide  and  a  set  of  planar  branch  waveg¬ 
uides,  one  per  chip,  each  of  which  couples  one  mode  out  of  the  trunk  waveguide.  The  elimination 
of  the  clock  skew  is  accomplished  by  taking  advantage  of  the  different  group  velocities  of  the 
modes  and  tailoring  the  propagation  constants  of  the  trunk  waveguide  according  to  the  location 
of  the  respective  chip  on  the  wafer. 

The  Darboux  transformations  of  the  inverse  scattering  theory  is  employed  to  design  refractive 
index  profiles  of  the  trunk  and  branch  waveguides,  using  the  set  of  propagation  constants  selected 
based  on  the  length  of  each  detector  from  the  source  and  the  group  velocity  of  the  mode  carrying 
the  dock  or  the  data  to  that  particular  detector  point.  Coupling  between  the  trunk  and  the  branch 
waveguides  are  also  analyzed.  It  is  theoretically  possible  to  ensure  nearly  100%  coupling  from  the 
trunk  to  the  branch,  although  the  trunk  and  the  branch  waveguides  are  not  identical.  Techniques 
for  ensuring  a  smooth  trunk  refractive  index  profile  are  investigated.  The  relationship  between 
circuit  size,  spread  of  the  propagation  constants,  and  allowable  circuit  loss  are  examined  in  detail. 
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3)  Analysis  of  Coupling  in  Multilayer  Planar  Optical  Waveguides 
using  Inverse  Scattering  Theory 

Coupling  between  waveguides  in  a  multilayer  structure  is  the  cornerstone  of  optical  spatial 
switching.  The  traditional  weak  ooupling  analysis  of  interacting  waveguides  has  been  reformulated 
in  the  language  of  scattering  theory.  We  show  that  the  coupling  coefficients  describing  the 
interaction  of  two  neighboring  waveguides  have  straight  forward  representations  in  terms  of  their 
scattering  data,  eliminating  the  need  to  explicitly  calculate  the  field  dependent  interaction  integrals, 
and  replacing  the  integrals  with  straightforward  algebraic  expressions  involving  the  guided  mode 
propagation  constant  and  the  residue  of  the  reflection  coefficient. 

(4)  Inverse  Scattering  Theory  for  the  Design  of  Optical  Waveguides 
with  Same  Propagation  Constants  for  Different  Frequencies 

We  have  developed  inverse  scattering  theory  for  designing  planar  optical  waveguides  pos¬ 
sessing  the  same  prescribed  propagation  constants  for  different  transmission  frequencies.  The 
design  problem  for  TE  modes  is  transformed  and  reformulated  to  an  equivalent  inverse  problem 
for  Schrodinger  equation.  Then  using  inverse  scattering  theory,  the  potential  as  a  function  of  a 
modified  spatial  variable  is  recovered.  Next  the  important  problem  of  finding  an  explicit  relation 
between  the  actual  spatial  variable  and  the  modified  spatial  variable  is  solved  and  a  systematic 
procedure  is  developed  for  designing  waveguides  which  have  the  same  propagation  constant 
for  different  light  frequencies.  Existence  and  uniqueness  questions  are  also  studied  and  design 
examples  are  worked  out. 

(4)  Design  of  Optical  Fibers  with  Same  Propagation  Constants  for  Different 
Azimuthal  Modes  and  Its  Application  to  Image  Transmission 

Inverse  scattering  theory  for  fixed  angular  momentum  has  been  modified  and  applied  to  the 
design  of  optical  fibers  possessing  prescribed  propagation  constants  for  different  azimuthal  modes. 
Such  designed  optical  fibers  find  application  in  the  transmission  of  spatial  images.  The  dephasing 
in  the  spatial  images  which  is  generally  attrfouted  to  the  different  modal  dispersion  of  different 
modes  can  be  totally  eliminated.  The  results  obtained  using  inverse  scattering  theory  has  been 
validated  using  a  finite  difference  base  direct  procedure.  Also,  sensitivity  studies  on  the  refractive 
index  profiles  have  been  investigated. 
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(5)  Finite  Difference  Analysis  of  Optical  Waveguides 

To  test  the  validity  of  any  design  obtained  using  the  inverse  scattering  theory,  an  alternative 
independent  approach  is  necessary.  So,  we  have  developed  a  finite  difference  based  direct 
scattering  solver  that  can  verify  the  design  results  obtained  using  the  inverse  scattering  theory. 

We  have  formulated  a  matrix  eigen  value  problem  for  cylindrical  and  planar  optical  waveguides 
from  a  set  of  finite  difference  equations.  Numerical  solution  of  this  problem  yields  the  propagation 
constants  for  propagating  modes.  The  method  can  be  used  for  arbitrary  index  profiles,  does  not 
require  the  explicit  evaluation  of  Bessel  or  modified  Bessel  functions,  and  does  not  use  iterative 
methods  to  search  for  the  propagation  constants  as  was  the  case  in  earlier  proposed  methods 
using  finite  differences.  The  method  we  have  developed  is  accurate,  fast  and  simple.  We  have 
established  the  convergence  and  stability  of  this  method,  and  explored  the  effects  of  finite  cladding 
width  on  the  dispersion  characteristics.  The  computer  software  we  have  developed  can  be  used 
interactively  to  explore  dispersion  in  optical  fibers  and  optical  waveguides. 

(6)  Nonlinear  Optical  Waveguides 

The  knowledge  gained  in  the  process  of  applying  spectral  inverse  scattering  theory  to  the 
design  optical  devices  is  useful  in  extending  to  the  study  of  nonlinear  optical  waveguides.  The 
problems  that  limit  the  information  carrying  capacity  of  the  soliton  based  communication  are  the 
soliton-soliton  interaction  and  the  soliton  self  frequency  shift  known  as  Gordon-Haus  effect.  We 
are  exploring  both  the  effects.  Also  we  are  studying  nonlinear  periodic  and  aperiodic  media  both 
from  the  inverse  scattering  point  of  view  and  using  simulation.  The  goal  is  to  synthesize  optical 
bistable  devices  which  can  be  extended  to  the  design  of  optical  logic  devices.  Optical  logic  devices 
are  crucial  to  the  development  of  *->ntical  computers. 


4 


LIST  OF  PUBLICATIONS  RELATED  TO  THE  PROJECT 


1 .  D.  W.  Mills  and  L.  S.  Tamil,  “Coupling  in  Multilayer  Optical  Waveguides:  An  Approach  Based 
on  Scattering  Data,”  submitted  to  J.  Light  Wave  Tech  not. 

2.  L.  S.  Tamil  and  Y.  Lin,  “  Synthesis  of  Optical  Ffoers  with  Same  propagation  Constant  for 
Different  Azimuthal  Modes:  Application  to  Imaging,”  under  preparation. 

3.  L.  S.  Tamil  and  Yun  Lin,  “  Synthesis  and  Analysis  of  Planar  Optical  Waveguides  with 
Prescribed  TM  Modes,”  J.  Opt.  Soc.  Am.  A,  in  print. 

4.  D.  W.  Mills  and  L.  S.  Tamil,  “  Synthesis  of  Guided  Wave  Optical  Interconnects,”  IEEE  J. 
Quant.  Electron.,  in  print 

5.  G.  H.  Aicklen  and  L.  S.  Tamil,  “  Interactive  Analysis  of  Propagation  in  Optical  Rbers,” 
Computer  Appl.  Engrg.  Edu.,  Vol.  11,  No.  3,  pp.  197-204  (1993). 

6.  L  S.  Tamil  and  G.  H.  Aicklen,  “Analysis  of  Optical  Fibers  with  Arbitrary  Refractive  Index 
profiles:  Accuracy,  Convergence,  and  Effects  of  Finite  Cladding ,*  Opt.  Commun.,  Vol.  99, 
pp.  393-404,  June  15,  1993. 

7.  D.  W.  Mills  and  L.  S.  Tamil,  “  A  New  Approach  to  the  Design  of  Graded-lndex  Guided  wave 
Devices,”  IEEE  Microwave  and  Guided  Wave  Lett.  Vol.  1,  No.  4,  April  1991. 

8.  D.  W.  Mills  and  L.  S.  Tamil,  “Analysis  of  Optical  Waveguides  Using  scattering  Data,”  J.  Opt. 
Soc.  Am.  A..  Vol.  9,  No.  10,  pp.  1769-1778,  Oct.  1992. 

9.  L.  S.  Tamil  and  Y.  Yu,  “  A  Beam  Propagation  Technique  to  Analyze  Integrated  Photonic 
Circuits,"  Microwave  and  Opt.  Techno).  Lett.,  Vol.  5,  No.  12,  pp.  617-621,  November  1992. 

10.  L.  Tamil,  “  Synthesis  of  Guided  Wave  Optical  Interconnects,”  Conf.  on  Emerging  Optoelec¬ 
tronic  Techno!.,  Bangalore,  India,  Dec.  16-20,  1991. 

11.  “Synthesis  of  Integrated  Optical  Devices:  An  Inverse  Scattering  Approach,  Progress  in 
Electromagnetics  Research  Symposium,  July  12-16,  1993,  Pasadena,  CA. 

12.  D.  W.  Mills,  *  Guided  Wave  Optical  Interconnects:  An  inverse  Scattering  Approach,”  Ph.D. 
Dissertation,  University  of  Texas  at  Dallas,  TX,  Aug.  1992. 

13.  Y.  Lin,  “  Synthesis  of  Planar  Optical  Waveguides  with  Preserved  TM  Modes,”  M.  S.  Thesis, 
University  of  Texas  at  Dallas,  TX,  Aug.  1992. 


5 


Appendix  A 


L  S.  Tamil  and  Y.  Lin,  “  Synthesis  and  Analysis  of 
Planar  Optical  Waveguides  With  Prescribed  TM  Modes” 

J.  Opt  Soc.  Am.  A.,  in  print 


6 


Synthesis  and  Analysis  of  Planar  Optical 
Waveguides  With  Prescribed  TM  Modes 


Lakshman  S.  Tamil  and  Yun  Lin 
Erik  Jonsson  School  of  Engineering  and  Computer  Science 

and 

Center  for  Applied  Optics 
University  of  Texas  at  Dallas,  Box  830688, 
Richardson,  TX  75083-0688 
(214)  690-2197 


7 


Abstract 


An  inverse  scattering  approach  to  designing  optical  waveguides  with  prescribed 
propagation  characteristics  of  TM  modes  is  presented.  The  refractive  index  profile  of 
the  waveguide  is  formulated  as  a  solution  to  a  nonlinear  differential  equation  whose 
forcing  function  is  the  potential  obtained  from  the  application  of  inverse  scattering  theory. 
This  method  can  reconstruct  smooth  refractive  index  profiles  for  planar  waveguides  that 
support  single  mode  or  multi-modes.  Both  the  cases  of  zero  and  non-zero  reflection 
coefficients  characterizing  the  transmission  properties  of  waveguides  are  discussed  here. 
A  direct  analysis  technique  based  on  finite  difference  scheme  has  been  formulated 
to  verify  the  results  obtained  by  inverse  scattering  method  and  they  are  in  excellent 
agreement 
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1.  Introduction 


The  conventional  method  of  designing  optical  waveguiding  structures  is  to  assume  a 
refractive  index  profile  and  solve  the  governing  differential  equation  to  find  the  various 
propagating  modes  and  their  propagation  characteristics.  If  the  propagation  characteristics 
do  not  meet  the  expected  behavior,  the  refractive  index  is  changed  and  the  propagation 
characteristics  are  again  evaluated;  this  is  repeated  until  the  expected  propagation  behavior 
of  the  modes  are  obtained.  This  being  an  iterative  procedure,  it  is  time  consuming.  Also, 
to  obtain  certain  arbitrary  transmission  characteristics,  one  may  not  be  able  to  imagine 
the  right  initial  refractive  index  profile.  It  is  important  to  understand  that  we  normally 
come  up  with  initial  profiles  that  have  mathematically  a  closed  form  such  as  parabolic, 
secant  hyperbolic  etc. 

The  procedure  discussed  in  this  paper  as  opposed  to  the  direct  method,  starts  with 
the  required  propagation  characteristics  of  the  waveguide  and  obtains  the  refractive  index 
profile  as  the  end  result  This  is  achieved  by  transforming  the  wave  equation  for  both 
the  TE  and  TM  modes  in  the  planar  waveguide  to  a  Schrodinger  type  equation  and  then 
applying  the  inverse  scattering  theory  as  formulated  by  Gelfand,  Levitan  and  Marchenko 
[1-2].  The  inverse  scattering  problem  encountered  here  has  a  direct  analogy  to  the 
inverse  scattering  problem  of  the  quantum  mechanics.  The  refractive  index  profile  of  file 
planar  waveguide  is  contained  in  the  potential  of  file  Schrodinger  type  equation  and  the 
propagating  modes  are  the  bound  states  of  the  quantum  mechanics  [3]. 

An  inverse  scattering  theory  with  zero  reflection  coefficient  characterizing  the  prop¬ 
agation  property  has  been  applied  earlier  to  the  design  of  planar  waveguides  by  Yukon 
and  Bendow  [43.  In  that  work,  the  refractive  index  profiles  were  constructed  only  for 
the  prescribed  TE  modes.  The  inverse  problem  of  designing  optical  waveguides  whose 
transmission  property  is  characterized  by  a  non-zero  reflection  coefficient  has  been  solved 
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for  TE  modes  by  Jordan  and  Lakshmanasamy  [5],  In  this  paper  we  have  applied  die  in¬ 
verse  scattering  theory  with  both  the  zero  and  the  non-zero  reflection  coefficient  to  design 
planar  waveguides  with  prescribed  TM  modes. 

In  section  2  we  review  the  problem  of  die  electromagnetic  wave  propagation  in  a 
planar  waveguide  for  both  TE  and  TM  cases  [6],  then  we  present  a  way  to  transform 
wave  equations  to  Schrodinger  type  equations.  In  section  3,  we  review  the  Kay’s  inverse 
scattering  theory  [7]  and  Gelfand,  Levitan  and  Marchenko  equation  [1-2].  The  inverse 
scattering  theory  is  then  applied  to  planar  waveguides  for  the  case  of  TM  modes  in  the 
zero  and  the  non-zero  reflection  coefficient  conditions  separately.  The  single  mode  and  the 
multi-mode  refractive  index  profiles  with  prescribed  TM  modes  are  obtained  by  solving  a 
nonlie ar  differential  equation  using  the  Runge-Kutta’s  fourth  order  approximation  method 
as  discussed  in  sections  4  and  5.  Construction  of  the  potentials  for  a  single  mode 
planar  waveguide  for  the  non-zero  reflection  coefficient  case  using  rational  function  of 
wavenumber  as  reflection  coefficients  is  presented  in  section  5. 

In  order  to  verify  the  results  obtained  by  inverse  scattering  theory  we  have  developed 
an  efficient  finite  difference  method  to  find  the  propagation  constants  of  guided  TE 
and  TM  modes  and  is  presented  in  section  6.  We  start  from  the  wave  equations  for 
TE  and  TM  modes  and  transform  them  to  a  set  of  finite  difference  equations.  Then 
a  matrix  eigenvalue  equation,  from  which  the  propagation  constants  can  be  found,  is 
constructed.  The  numerical  results  are  obtained  for  several  graded-index  waveguides, 
and  we  have  compared  these  results  to  the  previously  published  analytical  solutions  and 
results  obtained  by  other  numerical  methods.  The  conclusions  are  given  in  section  7. 
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2.  Physical  Model  of  a  Planar  Waveguide 

The  wave  equations  for  the  inhomogeneous  planar  optical  waveguides  can  be  derived 
from  the  Maxwell’s  equations.  If  we  take  z  as  the  propagation  direction  and  let  u> 
represent  the  frequency  of  laser  radiation,  we  have  the  following  wave  equations  for  one 
dimensional  inhomogeneous  planar  waveguides  [6] 

+  [«&(*)  -  £,(x)  =  0  (1) 

for  TE  modes  and 

for  TM  modes.  The  planar  waveguide  we  are  considering  here  has  a  refractive  index 
which  varies  continuously  in  the  x  direction.  For  the  planar  optical  waveguide  shown 
in  Fig.  1,  our  problem  is  to  find  the  refractive  index  profile  function  in  the  core  for  a 
set  of  prescribed  propagation  constants. 

We  assume  that  this  planar  waveguide  has  a  refractive  index  profile  guiding  N  modes. 
The  propagation  constants  {/?„}  are  korii  >  fii  >  fh  >  —0N  >  fco«oo»  in  which  rtoo 
is  the  value  of  n(x)  as  x  — ►  oo  and  ni  =  supn(x).  Designing  an  optical  waveguide  is 
analogous  to  the  inverse  problem  encountered  in  quantum  mechanics.  We  are  trying  to  get 
the  potential  function  from  the  given  bound  states  and  scattering  data.  The  wave  equation 
for  the  TE  modes  can  be  easily  transformed  to  an  equivalent  Schrodinger  equation 

^j£,(x)  +  [k2  -  V(x)]  £,(*)  =  0  (3) 

by  letting 

V(x)  =  -kl  [n2(x)  -  n^J  (4) 
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and 


k2  =  -4  =  -(Pi  -  4nL)  • 


We  can  see  in  our  case  the  potential  function  V(x)  is  continuous  and  V(x)  — »  0  as 
|  x  |— ►  oo.  The  TE  mode  cases  have  been  solved  by  Yukon  and  Bendow  [4]  and  Jordan 
and  Lakshmanasamy  [5],  so  our  discussion  will  be  restricted  to  TM  modes. 

We  now  need  to  transfer  the  wave  equation  for  the  TM  modes  to  Schrodinger  type 
equation  in  order  to  apply  the  inverse  scattering  method.  In  wave  equation  (2),  the  first 
derivative  of  Ex  can  be  eliminated  if  we  let  Ex(x)  —  e_1/2(x)$(x).  The  wave  equation 
then  becomes 


d$2  1  d2e(x) 

dx2  2e(x)  dx2 


3  (d4x)V 
4c2(x)  \  dx  J 


$  +  (A?oe(a:)  -  /92)$ 


=  0  . 


(6) 
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3.  Inverse  Scattering  Theory 

The  inverse  scattering  theory  of  Kay  and  Moses  [7]  provides  us  a  way  to  obtain  the 
potential  from  the  reflection  coefficient  which  characterizes  the  propagation  properties  of 
the  planar  waveguide.  As  the  potential  we  defined  vanishes  at  infinity,  we  can  apply  the 
Gelfand-Levitan-Marchenko  (GLM)  equation  to  solve  our  problem.  Let  us  consider  a 
time-dependent  formulation  of  the  scattering.  Taking  the  Fourier  transform  of  Eq.  (7); 
the  transform  pairs  are  $(x,fc)  ^(x,*)  and  k  t,  we  obtain 

-  V(x)*(x,t)  =  0  ,  (10) 

in  which  t  is  the  time  variable  with  the  velocity  of  light  c  =  1.  The  incident  plane  wave 
is  represented  by  the  unit  impulse 

\P(x,  t)  =  S(x  —  t)  ,  x  <  0  ,  t<  0,  (11) 

which  will  produce  the  reflected  transient  wave  function 

00  N 

R(x  +  t)  =  ±-  J  r Ane-'Ms+‘)  ,  (12) 

-00  n=1 

where  k2  =  —k2  are  the  discrete  eigenvalues  of  Schrodinger  type  equation.  (7),  r(k) 
is  the  complex  reflection  coefficient,  An  are  arbitrary  constants  normalizing  the  wave 
equation  such  that 

+oo 

J  $(x)$*(x)dx  =  1  .  (13) 

— OO 


The  reflected  transient  is  produced  only  after  the  incident  unit  impulse  has  interacted 
with  the  inhomogeneous  core  of  the  optical  waveguide  and  therefore 

R(x  +  t)  =  0  for  x+f  <  0  .  (14) 
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A  linear  transform  independent  of  k  can  now  relate  the  wave  amplitude  $(x,  t)  in  the 
core  region  with  the  wave  amplitude  <)  in  the  exterior  region 


*(M) 


,  *o(*,0 +/*(*,«') *<>({',<)<«’  I>0 

k^o(x,f)  x<0 


Here  the  exterior  field  is 


(15) 


*)  =  6(z  —  t)  +  R(x  + 1)  .  (16) 

From  physical  consideration,  since  ^(x,t)  is  a  rightward  moving  transient 

$(x,<)  =0  for  t  <  x  .  (17) 

Thus  the  kernel  K(x,t)  =  0,  for  t  >  x  and  K(x,t)  =  0  for  t  <  —x.  Substituting  Eq. 
(16)  into  Eq.  (15)  and  using  Eqs.  (14)  and  (17)  yield  the  integral  equation 

X 

K(x,t)  +  R(x  + 1)  +  J  k(x,^r(^  =0  t<x.  (18) 

—X 

We  can  show  that  by  substituting  Eq.  (15)  into  Eq.  (10)  the  kernel  K(x,t)  satisfies 
a  differential  equation  of  the  same  form  as  Eq.  (10)  provided  the  following  conditions 
are  imposed 


K(x,-x)=Q  ,  (19) 

and 

2±K(x,x)  =  V(x)  .  (20) 

We  now  could  see  how  the  solution  of  the  integral  Eq.  (18)  for  the  function  K(x,t)  can 
lead  to  the  synthesis  of  optical  waveguides. 
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4.  Design  Example  1:  Zero  Reflection  Coefficient 

The  reflection  coefficient  characterizes  the  propagation  properties  of  the  optical 
waveguides.  The  zero  reflection  coefficient  characterizes  a  system  with  propagating 
modes  only  where  as  the  non-zero  reflection  coefficient  characterizes  a  system  with  both 
guided  and  nonguided  modes.  Let  us  first  consider  the  special  case  of  zero  reflection 
coefficient  [8].  Substituting  Eq.  (12)  for  r(k)  =  0  in  GLM  equation  (18)  we  have 

N  N  » 

+  I  K(x, =  0  .  (21) 

n=l  n=l  4,, 

It  is  clear  from  the  above  equation  that  the  solution  for  K(x,  t)  should  have  the  form  [8] 

N 

*(*,!)  =  £/.(*)««-'  .  (22) 

n==  1 

Substituting  Eq.  (22)  into  Eq.  (21)  produces  a  system  of  equations  for  /„(x): 


An  S  ( *"+**'  j + fn{x 


)  +  AneKnX  =  0 


where  n  =  1, 2,  ...N,  This  system  can  be  conveniently  written  as 

[A]  [f ]  +  [B]  =  0  (24) 

where  [f]  and  [B]  are  column  vectors  with  /„,  and  Bn  =  An  exp  (nnx)  respectively,  and 


A  is  a  square  matrix  with  elements 


Avn  —  8vn  +  At, 


e(lt»+Kn)x ' 

Kv  +«n 


in  which  6„n  is  a  Kronecker  delta.  The  solution  for  f  is  f  =  —A  XB  and  then  from  Eq. 
(22)  K(x,x)  =  E^f  where  E  is  the  column  vectors  with  element  En  —  exp(/c„x)  and 


T  denotes  transpose.  Now, 


-A™  =  Ave^+K^x  =  BnEn 
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and  so 


K(x,x)  —  Enfn  —  —EnAv^Bv  —  (27) 

when  written  with  subscript  notation  and  the  summation  convention.  The  K(x,x)  given 
by  Eq.  (22)  can  be  recognized  in  the  form 

K(x,x)  =  =  ^ln(detA)  (28) 

and  therefore  the  potential  V(x)  according  to  Eq.  (20)  is 

V(z)  =  -2-^r  In  (det A)  .  (29) 

dxz 

Given  N  modes  with  desired  propagation  constants,  we  can  obtain  a  potential  function 
as  given  by  Eq.  (29).  Here  we  have  N  degrees  of  freedom  due  to  N  arbitrary  constants 

{A„|n  =  l,2..JV}. 

For  TE  modes  the  refractive  index  profiles  is  simply  given  by 

Ax)  =  A  -  oo) 

K0 

in  which  Jfco  is  the  free  space  wave  number.  Where  as  for  TM  modes,  obtaining  the 
refractive  index  profile  is  little  more  complicated  as  it  is  a  solution  to  a  nonlinear  differ¬ 
ential  equation  (8).  The  nonlinear  differential  equation  can  only  be  solved  numerically. 
Equation  (8)  is  first  transformed  to  a  convenient  form  by  setting  e(z)  =  we  then 
obtain, 

(3i) 

This  is  a  constant  coefficient  equation  which  yields  the  refractive  index  profile  \J e(z) 
provided  the  potential  V(x)  is  given. 
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To  demonstra^  some  practical  examples,  let  us  compute  the  refractive  index  profiles 


for  two  cases:  the  single  mode  case  and  the  N  mode  case. 

For  the  single  mode  case,  Eq.(23)  becomes 

A, c*'*  +  /,(*)  +  =  0  . 

Then,  the  potential  has  the  form 

V(  1  -  ~^iAie2ltlX 
*  (1  +  Aie2*1* /2k\)2  ’ 

where  A\  is  an  arbitrary  constant  and  note  that  k\  can  be  obtained  from 

«1  =  ~  ^Woo  • 


(32) 


(33) 


(34) 


For  a  desired  propagation  constant  ft,  we  can  get  a  set  of  refractive  index  profiles 
corresponding  to  different  arbitrary  choice  of  A\,  see  Fig.  2.  We  use  the  following 
data  relating  to  waveguide:  n(oo)  =  n,  =  2.177,  wavelength  A  =  0.8 /im  and 
ft  =  17.20  (/im)-1.  We  obtained  the  refractive  index  profiles  by  solving  Eq.  (31)  using 
the  potential  V(x)  obtained  from  Eq.  (33).  Runge-Kutta’s  fourth  order  approximation 
is  applied  in  solving  the  differential  equation  (31)  [9].  We  can  see  from  Fig.  2  that  the 
maximum  value  of  refractive  index  lies  on  the  positive  side  of  x  =  0  when  A\  <  2«i; 
on  the  negative  side  of  x  =  0  when  A\  >  2k\  and  at  x  =  0  when  A\  =  2k\. 

On  substituting  A\  =  2k\  into  Eq.  (33)  yields 

V(x)  =  —2K,\sech2Kix  .  (35) 


This  potential  is  everywhere  negative  and  goes  to  zero  as  x  goes  to  infinity.  Also  the 
potential  is  symmetric  about  its  minimum  point  We  can  truncate  the  potential  at  the 
point  where  the  potential  is  1%  of  its  maximum  value  to  find  the  width  of  the  core  d. 
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The  refractive  index  profile  corresponding  to  this  potential  is  shown  by  continuous  line 
in  Fig.  2. 

Similarly,  for  the  N  mode  case,  we  need  to  construct  the  potential  first  using  Eq.  (25) 
and  then  solve  the  nonlinear  differential  equation  (31)  for  the  refractive  index  profiles. 
For  a  set  of  prescribed  propagation  constants,  every  arbitrary  choice  of  normalization 
constants  will  produce  a  different  potential  and  a  corresponding  refractive  index  profile. 
In  order  to  construct  a  symmetric  refractive  index  profile  with  single  peak,  we  found  that 
the  normalization  constants  {An  |  n  =  1, 2...N)  must  satisfy  the  following  equation  [10] 

An  =  y/2 KnPn  ,  (36) 


where 


Pn  =  (-1) 


n— 1 


N 

n 

i/=l 


Ky  +  Kn 
Ky  &n 


n  =  1,2,  ...N 


(37) 


for  the  reflectionless  case.  Here  N  is  the  number  of  guided  modes  in  the  planar 
waveguide.  For  the  case  N  =  5,  using  sets  of  arbitrary  normalization  constants 
{An  |  n  =  1,2 ...N}  we  have  computed  die  refractive  index  profiles  and  these  are  shown 
in  Fig.  3.  The  symmetric  profile  obtained  using  the  condition  (36)  is  shown  by  continuous 
line  in  the  figure. 


5.  Design  Example  2:  Non-Zero  Reflection  Coefficient 

In  the  previous  section,  we  took  advantage  of  assuming  that  the  reflection  coefficient 
is  zero,  which  simplified  the  problem  a  lot  Now  we  are  going  to  solve  the  problem  with 
non-zero  reflection  coefficient  We  follow  here  the  work  of  Jordan  and  Lakshmanasamy 
[5]. 

We  take  the  rational  function  approximation  for  our  scattering  data.  We  represent  our 
reflection  coefficient  using  a  three-pole  rational  function  of  transverse  wave  number  k  [5], 
the  three  poles  are:  one  pole  on  the  upper  imaginary  axis  of  the  complex  k  plane,  which 
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represents  discrete  spectrum  of  function  R(x  + 1)  (see  equation  (12))  characterizing  the 
propagating  mode,  and  two  symmetric  poles  in  the  lower  half  of  the  k  plane,  which 
represent  the  continuous  spectrum  of  R(x  + 1)  characterizing  the  unguided  modes.  The 
three-pole  reflection  coefficient  can  be  written  as 


(k-k1)(k-k2)(k~k3) 


(38) 


where  r0  can  be  determined  by  the  normalization  condition  r(0)  =  —1,  this  condition 
ensure  total  reflection  at  k  =  0.  k\,  fa  have  following  forms:  k\  =  —c\  —  *C2  and 
fc2  =  ci  —  tC2.  The  third  pole  on  the  positive  imaginary  axis  is  fc3  =  ia. 


The  pole  positions  are  confined  to  certain  “allowed  regions”  that  are  determined  by 
the  law  of  conservation  of  energy,  which  can  be  represented  by  |r(fc)|2  <  1  for  all  real 
A:;  see  Fig.  3  in  Ref.  [5]  for  details. 

It  has  been  shown  that  the  reconstructed  potential  function  V(x)  has  following  form 


[5]: 


V(x)  =  2 


4»rM_ar(l)A-(t)j(  A(*)_) 


dx  '  7  y  '  dx 

in  which,  a  and  b  are  column  vectors,  and  are  given  by 


A  x(x)b  , 


(39) 


and 


where 


and 


ar(x)  =  [1  x  eVlS  e~mx  e*1  e ~’hZ] 
br  =  [0  0  0  0  0  -  a(cf  +  c2a)]  , 


m 


'ia2  +  ci-c;+i(a2-4ci)1/2(<.2+4c?),/2 


(40) 

(41) 


(42) 


(43) 
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Matrix  A(x)  is  given  by 

'0  10  0  0 
0  0  f(rji)  a(cf  +  <4)  0 

0  0  0  0  /(  t?2) 

1  -x  e~'1lX  e*11*  eT1*1* 

0  —1  —T)ie~mx  J?ie’hr 

.0  0  i)\e~*x  i)\t*x  T)%e-*x 

where  the  function  f(x)  is  constructed  by 

f(x)  =  x3  +  (2 c2  -  a)x2  +  [c2  +  cj  -  2ac2]x  -  a(cj  +  cf)  -  (45) 


0 

0 


a(c?  +  c3) 


trnx 

foe1*2* 


(44) 


So,  it  is  possible  to  construct  the  potential  from  the  three  poles  of  reflection  coefficient 
using  the  above  equations.  We  choose  here  two  examples.  In  example  1,  the  poles 
are  determined  by  the  following  parameters:  a  =  1.0,  c\  —  0.8,  and  C2  =  0.499; 
example  2  has  different  unguided  modes  characterized  by  c\  —  0.  05,  =  0.1 

and  the  same  propagating  mode  characterized  by  a  =  1.0.  Figure  4  shows  the  plots  of 
potential  functions  for  examples  1  and  2.  In  the  example  2,  we  see  that  the  potential 
is  everywhere  negative. 

Figure  5  shows  the  refractive  index  profiles  for  TM  mode  in  both  the  above  discussed 
examples  obtained  by  substituting  the  potentials  into  the  nonlinear  differential  equation 
(31)  and  solving  for  i/c(x).  We  notice  that  a  depressed  cladding  is  obtained  in  example 
1  and  we  also  see  that  the  profiles  we  found  here  resemble  the  profiles  we  normally  find 
in  practical  optical  waveguides  [11]. 
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6.  Verification  by  Analysis 

In  order  to  verify  the  results  obtained  by  inverse  scattering  theory,  a  finite  difference 
based  analysis  scheme  is  developed  here.  Using  this  method  we  find  the  propagation 
constants  of  guided  TM  modes  of  an  optical  waveguide  with  arbitrary  refractive  index 
profile.  Owing  to  its  simplicity  and  flexibility,  this  method  is  proved  to  be  very  effective. 

Now  we  consider  a  symmetric  planar  waveguide.  For  the  TM  modes  we  have  [6] 


with  the  Hy  component  obeying  the  wave  equation 


n 


2 


d_ 

dx 


n2{x)k$)Hv{x)  . 


(49) 


For  the  one  dimensional  graded  index  planar  waveguide,  the  refractive  index  is  a 
function  of  x,  and  the  wave  equation  can  be  transformed  to 


(50) 


If  Hv  and  its  derivative  are  single  valued,  finite  and  continuous  function  of  x,  we  have 


the  following  finite  difference  approximation  to  the  differentials 

Hi+1  -  Hi_! 


dH_ 

dx 


2  h 


(51) 


and 

<?H  ~  Hi+l  -  2 Hj  +  H%—\ 
dx 2  fi2 


(52) 
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in  which  we  have  used  H  instead  of  Hy  for  simplicity.  We  have  J/,_i  =  H(x  —  h); 
Hi  —  H(x);  and  Hi+ 1  =  H(x  +  h),  in  which  h  is  the  distance  between  the  grid  points 
and  i  is  the  index  of  the  grid  point  Now  we  obtain  the  following  equation  by  substituting 
the  finite  difference  approximation  of  the  first  and  the  second  derivative  of  H  into  the 
wave  equation  (50). 


1  drii 
riih  dx 


1  drti\ 
U{h  dx  ) 


Hi+ 1 


=  0 


(53) 


in  which  ri{  =  n(ih),  the  value  of  dni/dx  is  the  derivative  of  refractive  index  n  at  x  =  ih. 


We  have  chosen  three  grid  points  in  each  region:  the  substrate,  the  film  and  die  cover 
for  the  purpose  of  illustration  (see  Fig.  6).  For  the  case  considered  here  the  boundary 
conditions  are 


Ho  =0 


(54) 


and 

Ho  =  0  ,  (55) 

that  is,  the  field  vanishes  at  the  ends  of  the  cladding.  Absorbing  boundary  condition 
should  have  been  more  appropriate,  however  it  is  not  used  here. 

We  can  write  finite  difference  equation  at  every  grid  point  from  :  =  1  to  t  =  7.  We 
use  function  /(t)  to  represent  the  derivative  of  n(t)  which  is  obtained  again  by  a  finite 
difference  approximation  and  is  denoted  by 


/(•)  =  ^  ■  (56) 

Note  that  /  goes  to  zero  in  the  substrate  and  in  the  cover  region.  We  have  at  i  =  1, 
Hq  =  0  and  so 


(nM-/92-l)j/,+pff2  =  0  . 


(57) 
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At  i  =  2, 

^  ~  ^ -^2)  H2  +  ^Hi  =  0  .  (58) 

In  the  film,  we  have  at  t  =  3, 

-  ^/(1))*4 = °- (59) 

At  t  =  4, 

(h+iwm)H3+(n2{2)l4-fi  ~  T>)Ht+{*  ~  ^/(2))"5  -  «• (60) 

At  *  =  5, 

+  -  £)*+(£  ~  jffiS**))*  =  0.  (6D 

At  t  =  6, 

i«s  +  (njlg -/32-1)h!,  +  i//7  =  0  ,  (62) 

and  at  i  =  7,  since  He  =  0,  we  have 


i/f6  +  (»J*5  -  =  0 . 


(63) 


We  can  rewrite  these  finite  difference  equations  by  a  matrix  equation  for  conve¬ 
nience. 


‘an  -  P2 

°12 

0 

0 

0 

0 

0  ' 

H,' 

°21 

022  “  P2 

023 

0 

0 

0 

0 

H ; 

0 

«32 

033  “  P2 

034 

0 

0 

0 

Hz 

0 

0 

043 

a«-p2 

045 

0 

0 

H\ 

0 

0 

0 

054 

055  -  P2 

056 

0 

Hs 

0 

0 

0 

0 

065 

066  -  P2 

067 

H6 

0 

0 

0 

0 

0 

076 

077  -  p2. 

LflrJ 

(64) 
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in  which  the  elements  of  matrix  A  are  defined  by 


<*n  =  —^2  +  **2*0  =  aN,N  >  (65) 


<*12  : 

=  ^2  =  <**,*- 1  , 

(66) 

and  for  2  <  t  <  N  we  have 

2<i<Ni  ,Ni  +  N2<i<N 
Ni<i<Ni  +  N2 

(67) 

f  ^ 

<*.',. +i  =  l  _  /(») 

7?  n(i)h 

2<i<Ni  ,N1  +  N2<i<N 

Nx  <  i  <  Ni  +  N2 

(68) 

_  f  n2k0  ~  F 

M  "  {  n*(i)kl  -  £ 

2<t<iVi  ,  N\  +  N2  <  t  <  N 
Ni<i<Ni  +  N2 

(69) 

In  which,  f(i)  is  the  derivative  of  the  refractive  index  at  x  =  ih,  N\,  N2  and  N3  are  the 

number  of  grid  points  in  the  substrate,  film  and  cover  respectively,  and  N  =  Ni  +N2+N2 
is  the  total  number  of  grid  points.  The  other  elements  of  the  matrix  which  are  not  defined 
above  are  zeroes. 

The  matrix  A  can  now  be  split  into 

A  =  [B-/?2I],  (70) 

where  the  matrix  B  has  the  following  simple  form  and  I  is  the  identity  matrix. 


<*12 

0 

0 

<*22 

<*23 

0 

<*32 

<*33 

<*34 

0 

<*43 

<*44 

0 

0 

<*54 

0 

0 

0 

0 

0 

0 

0 

0 

0  ‘ 

0 

0 

0 

0 

0 

0 

<*45 

0 

0 

<*55 

<*56 

0 

<*65 

<*66 

<*67 

0 

<*76 

<*77  . 

(71) 
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The  matrix  equation  (64)  can  now  be  rewritten  in  the  form 


[B  —  /?2l]  H  =  0  . 


(72) 


To  find  the  propagation  constants  of  die  guided  TM  modes,  we  have  to  solve  the 
eigenvalue  problem  of  equation  (72),  which  has  a  nontrivial  solution  if  and  only  if  01 
are  eigenvalues  of  B.  So  finally  we  have 

{/))  =  VaiiBj  (73) 

for  both  odd  and  even  modes.  For  the  TE  modes,  the  situation  is  much  easier  since  the 
wave  equation  has  a  simpler  form  than  the  TM  modes.  The  field  component  Ew  obeys 
the  wave  equation  [6] 

-  n2(x)<3)  B,(x)  .  (74) 

We  can  find  a  matrix  expression  similar  to  the  one  we  found  for  TM  modes.  In  the  TE 
case  we  need  not  calculate  the  derivative  of  the  refractive  index  profile. 

For  the  given  refractive  index  profile  distribution  n  =  n(x)  the  matrix  B  can 
be  constructed  and  the  propagation  constants  {/?}  can  be  obtained  by  solving  for  the 
eigenvalues  of  this  matrix. 

Before  attempting  to  analyze  the  refractive  index  profiles  obtained  by  the  application 
of  the  inverse  scattering  theory,  we  want  to  see  whether  the  finite  difference  technique 
developed  here  provides  the  right  result  In  order  to  do  that  we  have  applied  the  technique 
to  various  refractive  index  profiles  such  as  parabolic  and  Gaussian  for  which  results  are 
already  available  in  the  literature  [12].  The  results  corresponding  to  TM  modes  arc  given 
in  Table  1  and  2.  They  all  show  that  our  analysis  technique  is  accurate  and  powerful. 
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Having  established  the  accuracy  of  the  finite  difference  technique,  now  we  can  use 
this  technique  on  the  arbitrary  refractive  index  profiles  we  have  obtained.  Figure  7 
shows  the  dispersion  characteristics  for  the  refractive  index  profile  with  single  symmetric 
peak  shown  in  Fig.  3.  The  normalized  frequency  V  has  been  determined  from  the 
waveguide  thickness  and  the  free  space  wavelength  of  the  propagating  modes,  here  we 
have  V  =  n\  —  n\  =  37.6883.  The  normalized  propagation  constant  we  used  here 
is  defined  by 


(75) 


Here,  n%  is  the  refractive  index  of  the  cladding  and  ni  is  the  maximum  refractive  index 
of  the  core.  As  we  see,  the  number  of  TM  modes  present  are  the  same  number  we 
started  with  in  reconstructing  the  profile.  The  refractive  index  profiles  corresponding  to 
the  non-zero  reflection  coefficient  as  shown  in  Fig.  5  when  analyzed  yields  a  dispersion 
characteristics  shown  in  Fig.  8.  Again  we  see  the  consistency  in  the  number  of  modes 
obtained  by  analysis  and  the  number  of  modes  used  in  the  synthesis  of  the  profile.  Though 
we  have  shown  that  the  number  of  modes  are  right,  it  is  not  a  sufficient  proof  that  the 
reconstructed  refractive  index  profiles  have  the  same  propagation  constants  for  each  of  the 
specified  modes.  To  check  this,  we  have  compared  (Table  3)  the  propagation  constants 
of  various  modes  we  used  in  reconstructing  the  refractive  index  profile  of  the  waveguide 
against  the  propagation  constants  obtained  by  analysis  for  the  normalized  frequency  at 
which  the  propagation  constants  are  prescribed.  We  see  that  last  two  columns  of  the 
table  agree  very  well.  This  shows  that  the  inverse  technique  outlined  here  can  be  used 
to  synthesize  waveguides  with  prescribed  TM  modes. 


26 


7.  Discussions  and  Conclusions 


We  have  developed  a  method  based  on  inverse  scattering  theory  that  can  be  used  to 
design  planar  optical  waveguides  which  transmit  prescribed  number  of  TM  modes  with 
prescribed  propagation  constants.  The  results  have  been  verified  using  finite  difference 
analysis.  This  procedure  in  conjunction  with  the  technique  for  designing  planar  optical 
waveguides  for  prescribed  TE  modes  developed  in  references  [4]  and  [5]  provide  the 
complete  inverse  scattering  procedure  for  designing  planar  optical  waveguides  with 
prescribed  propagation  characteristics.  However,  it  should  be  mentioned  that  only  the 
characteristics  of  one  of  the  kinds  of  the  modes  (IE  or  TM)  can  be  prescribed  in  a 
waveguide  as  they  are  governed  by  two  different  differential  equations. 

One  important  question  that  should  be  answered  when  we  fabricate  actual  waveguides 
with  refractive  index  profiles  obtained  using  the  technique  described  here  is  with  what 
precision  the  n(x)  should  be  fabricated  in  order  to  provide  the  desired  mode  configuration. 
To  answer  this  question  we  have  changed  V(x)  (V(x)  is  related  to  n(x)  through  Eq. 
(30)  )  uniformly  over  the  spatial  distance  x  by  1%,  5%  and  10%  and  have  computed 
the  corresponding  change  in  the  propagation  constants  for  a  typical  single  mode  profile. 
Figure  9  shows  the  plot  of  refractive  index  profile  corresponding  to  a  uniform  change  in 
V(z)  along  x  and  Table  IV  provides  the  computed  results  of  changes  in  the  propagation 
constant  due  to  the  uniform  change  in  V{x)  along  x.  We  see  that  a  change  in  AV/V 
in  the  range  of  1  to  5%  does  not  affect  significantly  the  mode  characteristics  of  the 
waveguide. 

It  is  also  important  to  analyze  the  effect  of  the  arbitrary  constants  An  on  the  shape 
of  the  resultant  refractive  index  profiles.  Figure  10  shows  the  variation  in  the  shape  of 
the  refractive  index  due  to  changes  in  the  choice  of  the  constants  {An  |  n  =  1,2,  ...N). 
Our  inference  is  that  the  shape  is  not  very  sensitive  to  die  changes  in  the  constants  An. 
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The  technique  developed  here  may  find  application  in  designing  waveguiding  struc¬ 


tures  for  spatial  transmission  of  images  and  optical  interconnections. 
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Table  1 


mode 

number 

7 

(3^/ko  from 
Ref.  [12] 

by  present  method 

N2=84 

N2=168 

0 

1.5966 

1.5966 

1.5966 

1 

1.5915 

1.5915 

1.5915 

2 

1.5864 

1.5864 

1.5864 

3 

1.5813 

1.5813 

1.5813 

4 

1.5762 

1.5762 

1.5762 

5 

1.5711 

1.5711 

1.5711 

6 

1.5659 

1.5658 

1.5659 

7 

1.5607 

1.5604 

1.5605 

8 

1.5556 

1.5546 

1.5548 

1.5503 

1.5498 

1.5499 

10 

1.5451 

1.5439 

1.5443 

11 

1.5399 

1.5387 

1.5390 

12 

1.5346 

1.5326 

1.5336 

13 

1.5294 

1.5288 

1.5290 

14 

1.5241 

1.5228 

1.5231 

15 

1.5188 

1.5139 

1.5143 
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Table  2 


mode 

number 

7 

p7/k0  from 
Ref.  [12] 

/37/ko  by  present  method 

N2=84 

N2=168 

0 

1.5984 

1.5984 

1.5984 

1 

1.5925 

1.5925 

1.5925 

2 

1.5867 

1.5868 

1.5867 

3 

1.5811 

1.5812 

1.5811 

4 

1.5756 

1.5757 

1.5757 

5 

1.5702 

1.5704 

6 

1.5649 

1.5651 

1.5650 

7 

1.5596 

1.5598 

1.5596 

8 

1.5545 

1.5542 

1.5544 

9 

1.5494 

1.5487 

1.5490 

10 

1.5444 

1.5424 

1.5431 

11 

1.5390 

12 

1.5341 

13 

1.5297 

1.5283 

1.5284 

14 

1.5251 

1.5220 

1.5222 

15 

1.5204 

1.5192 

1.5195 
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Table  3 


number  of 

mode 

modes 

number  7 

N=1 

0 

0 

N=2 

1 

0 

N=3 

1 

2 

0 

1 

N=5 

2 

3 

4 

0 

1 

2 

N=7 

3 

4 

5 

6 

prescribed  mode  spectra 

0y/ko  obtained  by  our 

Pylh 

analysis 

2.18997 

2.18995 

2.20556 

2.20553 

2.18417 

2.18398 

2.20926 

2.20916 

2.19140 

2.19100 

2.18061 

2.18036 

2.21288 

2.21266 

2.20003 

2.19968 

2.18998 

2.18968 

2.18278 

2.18254 

2.17845 

2.17797 

2.21466 

2.21452 

2.20473 

2.20449 

2.19630 

2.19606 

2.18927 

2.18915 

2.18397 

2.18379 

2.18010 

2.17997 

2.17778 

2.17753 
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Table  4 


prescribed 

effective 

index 

(0/ko) 

effective  index  (fi/ka)  obtained  by  our  analysis 

2.18997 

AV(x)  _  n% 

^1=5% 

=  10% 

2.18995 

2.18998 

2.19016 

2.19024 
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Synthesis  of  Guided  Wave  Optical  Interconnects 

Duncan  W.  Mills  and  Lakshman  S.  Tamil 

Abstract 

We  have  designed  a  guided  wave  optical  interconnect  which  reduces  or  eliminates  clock 
skew  by  ensuring  simultaneous  delivery  of  clock  pulses  to  chips  mounted  on  a  wafer.  The 
interconnect  consists  of  a  multimode  trunk  waveguide  and  a  set  of  branch  waveguides,  one  per 
chip,  each  of  which  couples  one  mode  out  of  the  trunk  waveguide.  The  elimination  of  dock 
skew  is  accomplished  by  taking  advantage  of  the  different  group  velodties  of  the  modes  inherent 
in  multimode  waveguides  and  suitably  tailoring  the  propagation  constants  of  the  trunk  waveguide 
according  to  the  location  of  the  respective  chip  on  the  wafer.  Inverse  scattering  theory,  specifically 
the  method  of  Darboux  transformations,  is  employed  to  design  the  refractive  index  profiles  of 
the  trunk  waveguides,  using  the  set  of  propagation  constants  selected  during  the  first  stag*  of  the 
design,  as  input  data.  It  is  shown  that  by  using  transverse  coupling  and  giir»Mr  design  of  the 
trunk  and  branch  waveguides,  efficient  coupling  from  the  trunk  to  the  branch  waveguides  can  be 
ensured.  Techniques  for  ensuring  a  symmetric  trunk  refractive  index  profile  arc  investigated. 
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1.  INTRODUCTION 


High-speed  computer  circuitry  requires  the  distribution  of  information  and/or  clock  pulses 
between  various  hardware  elements  within  the  system,  including  boards,  chips,  and  logical 
elements  within  a  chip.  Ideally,  the  clock  signals  reach  their  intended  destinations  simultaneously, 
but  in  practice,  the  exact  arrival  times  are  skewed  since  the  clock  signal  emanates  from  a  single 
source  to  various  locations  distributed  at  different  lengths  from  the  clock  source.  In  the  past,  a 
number  of  techniques  for  reducing  clock  skew  in  standard  VLSI  systems  consisting  of  metal  or 
polysilicon  interconnects  have  been  suggested.  These  include  layouts  composed  of  equal-length 
lines[l],  or  breaking  the  chip  into  a  blocks,  each  with  an  internally  generated  high-frequency 
dock  controlled  by  a  low-frequency  chipwide  clock[2].  Aside  from  the  fact  that  it  is  not  always 
practical  to  arrange  circuit  elements  to  meet  these  physical  requirements,  a  large  amount  of  metal 
wiring  is  required  to  implement  these  schemes.  The  trend  towards  higher  data  rates,  resulting 
in  skew  which  is  a  larger  percentage  of  the  clock  pulse  duration,  has  exacerbated  the  situation. 
This  paper  presents  a  method  for  designing  guided-wave  optical  interconnects  with  reduced  dock 
skew,  applicable  in  a  chip-to-chip  or  intrachip  situation. 

The  potential  advantages  offered  by  optical  interconnections  over  standard  wire  or  polysilicon 
lines  are  discussed  in  a  good  review  article[3],  which  suggests  that  optics  can  alleviate  problems 
stemming  from  resistive  and  capacitative  loading  in  wire/poly  lines,  which  is  deleterious  to  the 
signal  amplitude  and  shape,  particularly  at  higher  frequencies. 

In  this  paper,  it  is  proposed  that  graded-index  guided  wave  interconnects  can  effectively  re¬ 
duce  clock  skew  by  suitable  design  of  the  refractive  index  profile[4].  This  design  is  accomplished 
by  properly  tailoring  the  propagation  constants  of  the  guided  modes  to  provide  equal  propagation 
times  to  a  set  of  detectors.  The  scheme  presented  in  tbis^paper  employs  several  optical  channels, 
each  having  a  different  refractive  index  profile.  This  includes  a  main  multimode  channel  and 
several  single-mode  waveguides  coupled  to  the  main  line.  Total  system  design  takes  into  account 
the  problem  of  clock  skew  as  well  as  efficient  coupling  between  the  trunk  and  branch  waveguides. 

Section  2  of  this  paper  describes  the  relation  between  clock  skew  and  the  guided-mode 
spectrum,  followed  by  a  description  of  the  proposed  optical  interconnection  layout.  Section  3 
provides  the  physical  model  of  the  optical  waveguide  used  as  lire  building  block  of  the  optical 


interconnect  circuitry.  The  refractive  index  profile  of  the  multimode  guide  is  then  carried  out 
using  an  efficient  reconstruction  algorithm  which  generates  a  refractive  index  profile  based  on 
the  guided  mode  spectrum  and  desired  coupling  characteristics  between  the  main  waveguide  and 
the  single-mode  guides.  These  are  discussed  in  sections  4  and  5,  respectively.  Design  examples 
are  provided  in  section  6,  leading  to  conclusions  in  section  7. 


2.  OPTICAL  INTERCONNECT  CIRCUIT 


The  interconnect  network  is  to  be  mounted  on  a  GaAs  wafer  (10.16  cm.)  in  diameter,  as 
shown  in  Fig.l.  The  goal  of  the  interconnect  circuit  is  to  deliver  a  pulse  from  the  source  to  each 
of  the  detector  points  on  the  wafer  simultaneously.  The  circuit  consists  of  N  detectors  or  chips 
at  points  connected  by  a  network  of  integrated  optical  waveguides  consisting  of  a 

N- mode  trunk  line  feeding  N  branches,  which  are  generally  of  a  different  design  from  the  trunk. 
A  pulse  impressed  upon  a  given  mode  travels  at  the  group  velocity  vg,  where 

vg  du~  p  c1' 

so  that  a  pulse  traverses  a  given  length  L  in 


T  = 


w  n£ 
P  c2 


sec. 


(2) 


Consider  a  clock  pulse  launched  into  the  interconnect  at  point  S.  The  time  for  a  given  mode  to 
propagate  from  5  to  a  designated  jP(m)  is 

Tm  =  ir~k  L(m)  5ec-’  rn=l,  (3) 

Pm  C  — 

where  L(m>  is  the  distance  from  the  source  to  P(m),  and  pm  is  the  propagation  constant  of  the 
mode  delivering  the  signal  to  P(m)- 

For  the  purposes  of  this  analysis,  the  important  quantities  are  the  total  distances  from  the 
source  to  the  points  F(m).  Arranging  these  in  order  of  increasing  length. 


fyv  >  Ln-i  >  >  L\t 


(4) 


5! 


so  that  Ln  =  max  {L^},Li  =  min  {l/(m)},  the  points  P^)....P^  will  be  synchronized  if 
the  propagation  constants  satisfy 

02  _  ^2  03  _  ^3  0N  _  Ln 

fit  ~  Lt  ;  02  -  L2 ;  '"'0N-1  ~  Ijv-i'  1 

This  provides  us  with  a  set  of  JV  propagation  constants 

0t  <  02  <  ...  <  0N  (6) 

which  must  be  supported  as  propagating  modes  within  die  interconnect  The  propagation  constants 
themselves  are  restricted  to  the  range 

/joTl2  ^  0m  ^  koHmaxi  00 

where  fimax  is  the  maximum  core  refractive  index,  and  U2  die  refractive  index  of  the  cladding. 

The  design  of  the  interconnect  circuit  consists  of  two  interrelated  parts.  The  first  concerns 
the  design  of  the  refractive  index  profile  for  die  multimode  trunk  waveguide,  based  upon  the 
spectrum  generated  in  Eq.(5).  The  second  involves  die  design  of  the  branch  waveguides,  cadi 
of  which  must  efficiendy  couple  off  die  appropriate  mode  from  die  trunk  and  deliver  it  to  the 
designated  point.  This  raises  die  issue  of  waveguide  coupling.  In  sections  3  -  6,  we  illustrate 
the  application  of  inverse  scattering  theory  to  the  related  problems  of  refractive  index  profile 
design  and  coupling. 

3.  WAVEGUIDE  MODEL  FOR  OPTICAL  INTERCONNECT 

The  waveguide  model  consists  of  a  one-dimensional  planar  structure  with  graded-index  core 
n(z)  and  cladding  layers  of  constant  refractive  index,  as  shown  in  Ftg.2[5]. 

With  propagation  taken  along  the  z  axis,  the  TE  modes  take  the  form 

Ey(x)  (8) 

where  the  electric  field  Ey(x)  is  given  by 

-/?]£,  =  0.  (9) 
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Rearrangement  of  the  parameters  by  defining  the  transverse  wavenumber  k2  and  the  potential 


«(*)• 

k  =  klnf  -  p2 
v(x)  =  kUn2  -  n2(x)], 

brings  Eq.(2)  into  the  Schrodinger  equation  form. 


(10) 


<PEy 

dx2 


+  (k2  -  v(x))Ev  =  0. 


(11) 


Here,  ko  =  u/c  is  die  free  space  wavenumber,  /?  is  the  propagation  constant,  and  c  is  the 
velocity  of  light  in  vacuum.  From  these  considerations  it  is  dear  that  %  represents  the  asymptotic 
refractive  index  of  die  corresponding  waveguide,  provided 


v(x)  — ►  0  as  |x|  — *■  ±oo. 


(12) 


The  exact  model  for  die  waveguide  is  a  channel  geometry.  However,  for  the  sake  of  mathematical 
simplicity,  we  consider  the  planar  geometry  with  one  transverse  coordinate.  For  certain  separable 
refractive  index  profiles,  die  two-dimensional  refractive  index  surface  can  be  written  in  the 
additive  form[6]. 


«2(*>y)  =  n|  +  n2(x)  +  n2(y),  (13) 

in  which  case  the  y-depcndcnt  portion  of  the  refractive  index  can  be  designed  using  the  results 
of  planar  geometry,  resulting  in  a  complete  design  of  n(x,y). 

4.  RECONSTRUCTION  BY  TRANSFORMS 

Inverse  scattering  theory  provides  a  framework  whereby  die  potential  of  the  Schrodinger 

Mf  ?  ' 

equation  can  be  reconstructed  from  a  set  of  eigenvalues  selected  a  priori.  Inverse  reconstruction, 
based  on  the  solution  of  the  Gelfand-Levitan  integral  equation,  has  recently  been  applied  to  the 
design  of  monomode  waveguides[5],(7J.  In  general,  this  technique  is  cumbersome  when  several 
bound  states  are  present.  As  an  alternative,  we  will  employ  the  method  of  transformations  (known 
variously  as  Darboux  or  Cnim-Krein  transformations[8],[9])  to  obtain  multimode  potentials 
suitable  for  refractive  index  design  in  optical  interconnects. 


For  these  purposes,  it  is  useful  to  have  a  basic  understanding  of  scattering  parameters  related 
to  the  potentials  of  the  Schrodinger  equation,  which  we  assume  to  behave  asymptotically  as, 

v(x)  — ►  0,  x  — ►  ±oo.  (14) 

A  plane  wave  e+<**  incident  on  the  potential  from  x  —  -oo,  will  give  rise  to  a  reflected  portion 
taking  the  form, 

r_(  k)e~ikx  (15) 

as  i  -+  — oo,  as  well  as  a  transmitted  wave, 

t~(k)e+ikx  (16) 

as  x  — *•  oo  [10].  Similarly,  the  coefficients  r+(k)  and  t+(k)  can  be  defined,  where  it  can  be 
shown  that  t+(k)  =  t-(k)  =  t(k).  The  Schrodinger  equation  then  admits  a  pair  of  Jost  solutions, 
f+(k,x)  and  f-(k,x),  defined  according  to  their  asymptotic  behavior 

f±(k,x)e*,kx  =  1  x  — ►  ±oo,  (17) 

The  pairs  {/+(£,  *),/+(-fc,x)}  and  {/_(£, x),  f-{-k,  x)}  comprise  sets  of  linearly  independent 
solutions  to  the  Schrodinger  equation,  allowing  construction  of  the  linear  combinations  in  terms 
of  the  transmission  coefficient  t(k)  and  die  pair  of  reflection  coefficients  r+(fc)  and  r_(fc): 

/+(*,*)  =  /-(-*,  *)+^ /-(*•*),  (18) 
and 

/-(M)  =  ^A(-*,*)  +  ^/+(M).  09) 

The  scattering  data,  consisting  of  these  scattering  coefficients  and  the  bound  state  eigen¬ 

values, 

km  =  it cm  («m  >0),  (20) 

along  with  the  normalization  constants,  completely  characterize  the  form  of  the  potential.  The 
bound  state  wavefunctions  are  characterized  by  exponential  decay  for  large  |*|,  and  there  is  a 


54 


direct  one-to-one  correspondence  between  these  bound  states  and  the  guided  waveguide  modes 
characterized  by  a  discrete  spectrum  of  propagation  constants 


Pm  =  \Jk%nf  -  k£  =  \Jklnf  +  K.I,.  (21) 

It  is  clear  from  Eqs.  (17)-(19)  that  the  eigenvalues  are  poles  of  the  transmission  coefficient  which 
lie  on  the  upper  imaginary  axis  of  the  complex  k  plane  and  that  the  bound  state  wavefunctions, 
which  behave  asymptotically  as  are  merely  the  Jost  solutions  evaluated  at  these  poles: 

/±(*«m,*)  =  yf‘- F  (**»»  >  *  )>  (22) 

from  which  follows  the  important  relation, 

r+(*Km) 

A  rather  extensive  derivation  leads  to  alternative  representations  of  the  ratio  of  scattering 
coefficients  implied  by  Eq.(22),  useful  for  quantifying  waveguide  coupling[ll]: 

+oo 

2ik  ^  =  I  U(k,x)  v(x)  e*ik*  dx.  (24) 


Equation(5)  provides  a  prescription  for  constructing  a  spectrum  beginning  with  the  highest 
mode.  Pi.  whose  value  is  arbitrary,  subject  only  to  the  requirement  Pi  >  /^n2,  and  building 
upon  it  until  the  fundamental  mode,  characterized  by  Ph,  is  added  to  the  spectrum. 

The  transform  procedure  is  a  technique  which  allows  for  the  construction  of  N  -  mode  poten¬ 
tials  by  specifying  a  priori  a  set  of  bound  state  eigenvalues,  derived  from  the  set  {Pi,P?, Pn) 

~r  ■? 

via  Eq.(10): 


km  €  (25) 

where 

*N  >  kn-i  >  •••  >  «i  >  0.  (26) 
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The  following  discussion  does  not  constitute  proof  of  the  transformation  method,  for  which 
interested  readers  are  referred  to  reference  [9].  Rather,  it  oudines  the  practical  steps  necessary  to 
construct  a  particular  class  of  potentials  which  suit  the  purposes  of  refractive  index  profile  design. 

In  the  method  outlined  here,  we  are  adding  the  N  bound  states  to  some  chosen  initial 
potential  designated  t^(x)  which  is  assumed  to  contain  no  bound  states.  As  additional  input  data, 
we  require  die  explicit  form  of  the  Jost  solutions  f±{k,x)  associated  with  vo(ar). 

Defining  the  N  linear  combinations  of  the  Jost  solutions  as 

7m  =  (~l)m+1f+(iKm,x)  +  Pmfi(iKm,x)  T71  =  1, AT,  (27) 

where  pm  is  an  arbitrary  positive  definite  parameter,  the  corresponding  N-  mode  potential  is 
simply: 

eP 

»*(*)  =  Mx)  (28) 

In  the  above  equation,  the  quantities  W()  denote  Wronskians,  i.e., 

7l  72  .  .  IN 

7x  72  •  •  Yn 

W^(7i>72>—>7n)  =  •  .  .  .  .  ,  (29) 

AN-l)  (N- 1)  ■  '  (N- 1) 

7i  72  •  •  7m 

whose  rows  consist  of  functions  71. ..7 n  differentiated  with  respect  to  x  from  zero  to  TV  —  1 
times.  Equation(29)  clearly  illustrates  the  manner  in  which  the  N  bound  states  are  progressively 
added  in  stages  represented  by  rows  and  columns  of  the  Wronskian. 

The  Jost  solutions  corresponding  to  the  generated  potential  take  the  form 


(31) 


,,,,  k  -f  tKi  k  +  iK2  k  +  itiN 

r±(A)  =  (-1)"  i±£ 1  i±is....i±4^  «*<*), 

k  —  tKi  k  -  tK2  k  -  t Kpf 

cleaiiy  illustrating  the  presence  of  N  poles  representing  the  N  bound  states. 
The  normalization  constants, 


/£(fcm,a:)<i*  =  ImiResr+tkn)}, 


dm=\  fl(km,x)dx\  =  Im{Res  r_(fcm)}, 


which  play  a  critical  role  in  waveguide  coupling,  are  related  through  the  transmission  coefficient: 

cmdm=  (33) 

The  normalization  constants  also  transform  in  a  controlled  way  as  the  potential  is  constructed. 
It  can  be  shown  that  (see[9])t 


2  _  **  '-TO  p 

cm  —  .  rm  i 

Pm 


where 


\  i=  m  / 


m  =  l,2,..,iV. 


In  effect,  to  every  set  of  JV  3-tuples 


{®o(*)i*m!?m}  til  tx_1j  •••»  N 


there  corresponds  a  unique  N-  mode  potential  t?Ar(xj.  In  this  paper,  we  will  assume 

vo(x)  =  0, 


with  pertinent  scattering  data 


/*(*,x)  =  «**', 


m  =  i, 
R±(k)  =  0. 
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To  demonstrate  the  procedure  used  here,  we  have  used  three  sets  of  data  shown  in  Table  1, 
and  have  reconstructed  the  refractive  index  profiles  shown  in  Figs.(3)-(5).  Unless  otherwise 
indicated,  we  assume 


n  2  =  3.0 
A  =  0.9  nm , 


(39) 


throughout  this  paper.  In  Fig.3,  the  eigenvalues  arc  equally  spaced  and  all  pm  =  1.  The  resulting 
refractive  index  profile  is  symmetric  about  the  origin  and  provides  a  smooth  single  channel. 
If  near-degeneracies  are  introduced  into  the  spectrum,  splits  will  occur  in  the  refractive  index 
profile.  The  nature  and  extent  of  the  splits  will  depend  upon  the  number  of  modes  involved. 
This  is  illustrated  in  Fig.4,  where  we  introduce  a  quasi-degeneracy  across  the  entire  spectrum, 
causing  the  expected  split  of  the  refractive  index  profile  into  five  channels.  The  profile  remains 
symmetric  about  the  origin.  When  the  original  spectrum  is  restored,  but  one  or  more  of  the 
{ pm|m  =  1,2,  ...jiV}  deviate  from  unity  (Fig.(5)),  the  result  is  a  splitting  despite  the  wide 
spacing  of  the  eigenvalues. 

This  example  illustrates  the  application  of  the  reconstruction  procedure  to  typical  spectra 
compatible  with  GaAs  technology,  with  resulting  refractive  index  profiles  which  are  symmetric 
and  well-behaved.  It  is  clear  from  our  analysis  that  a  material  such  as  AlGaAs  [12],  with 
a  large  variation  in  refractive  index  as  a  function  of  mole  fraction,  is  well  suited  to  the 
proposed  interconnect  since  it  allows  for  a  relatively  laige  spread  of  propagation  constants,  and 
consequently,  greater  flexibility  in  placement  of  chips  on  the  wafer. 

As  a  check  of  our  algorithm,  one  hundred  data  points  representing  the  value  of  the  refractive 
index  profile  shown  in  Fig.3  were  used  as  input  to  a  finite-difference  algorithm,  solving  Eq.(9), 
whose  output  consisted  of  the  corresponding  five  propagation  constants.  Results  of  this  direct 
solution  are  given  in  the  second  column  of  Table  2,  showing  excellent  agreement  with  the  exact 
propagation  constants. 

Completion  of  the  optical  interconnect  problem  involves  design  of  the  branch  waveguides. 
In  the  next  section  we  discuss  transverse  coupling  between  the  trunk  and  arms,  and  show  how 
the  parameter  pm  can  be  adjusted  to  provide  the  desired  coupling  characteristics. 
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5.  COUPLING  TO  BRANCHES 

Efficient  power  transfer  from  the  trunk  to  each  branch  can  be  accomplished  by  transverse 
coupling  which  will  occur  over  the  interaction  length  represented  in  Fig.l  by  the  short  spans 
along  which  the  branches  are  parallel  to  the  trunk  waveguide.  The  analytical  approach  used  here 
is  based  upon  die  standard  perturbation  technique,  under  an  assumption  of  weak  coupling[13]. 
The  novelty  of  our  analysis  lies  in  relating  the  scattering  analysis  of  the  previous  section  to  the 
calculation  of  transverse  coupling  coefficients  derived  from  coupled  mode  theory. 

Figure  6  shows  two  neighboring  (non-overiapping)  waveguides,  each  of  which  is  assumed  to 
have  a  graded-index  core. with  refractive  index  profiles  n&(z)  and  nR(x),  respectively.  We  will 
assume  that  each  waveguide  separately  supports  y—  polarized  TE  modal  fields  El(z)  and  Er(x) 
with  propagation  constants  Pl  and  respectively.  The  interaction  between  the  two  fields  will 
be  represented  by  a  z—  dependent  linear  combination  of  the  individual  waveguide  modes: 

S(x,  z,  t)  =  A(z)  ER(x)e^-^  +  B(z)  EL(x)e^-0^ ,  (40) 


where  the  exact  form  of  the  z-  dependent  weighting  coefficients  A(z)  and  B(z)  are  to  be 
determined.  We  will  assume  that  the  coupling  is  weak,  i.c.. 


d?A(z) 


dz 2 


<< 


Pr 


dA(z) 


dz 


d*B(z) 


dz* 


<< 


0L 


dB(z) 


dz 


(41) 


The  interaction  between  the  waveguides  is  represented  by  first-order  differential  equations: 


dA(z) 

dz 


iKRLB(z)€~i^-P^X, 


(42) 


and 


The  coupling  coefficients. 


and 


=  A(z)S0t-p*>. 


KRL~  {2^^}’ 


(43) 


(44) 


(45) 
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are  defined  in  terms  of  the  various  integrals 


Tw 

Nrr  =  J  ER(x)  dx, 

-oo 

-foo 

Irl  =  J  Er(x)vr(x)Ei(x)  dx, 


TW 

iVL£,=  /  ££(x)dz 

£ 

Ilr=  J  EL(x)vL(x)ER(x)dx 

— OO 

The  coupling  coefficients  may  be  written  in  terms  of  the  scattering  data  as  follows.  Let  the  model 
consist  of  two  waveguides  described  by  potentials  Vf,(x)  and  vR(x).  For  the  purposes  of  the 
analysis,  let  us  shift  v^(x)  in  the  negative  x  direction  by  an  amount  s  so  that 

Er(x )  =  /*(*«>*) 

El(z)  =  /+(*'*,*  +  *) 

(48) 

vr(x)  =  %  H  “  nfi(x)] 
vl(x  +  5)  =  kl  [n^  -nj,(x)]. 

This  form  of  the  fields  was  chosen  so  that  they  have  the  simple  asymptotic  forms 

ER{x)  — ►  eKX  as  x  — ♦  -oo, 

(49) 

El{x)  -  as  z  —  oo. 

The  separation  s  is  arbitrary,  subject  only  to  the  condition  that  the  potentials  do  not  overlap  to 
any  large  extent. 

In  the  region  comprising  vR(x),  El(x),  the  field  of  the  lefthand  potential  taken  alone,  takes 
the  simple  form, 

El(x)  =  ,  (50) 
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enabling  the  interaction  integral  to  be  written  as 


= 


+OO 

J  f-(k,  x)  vr(x)  etkx  dx |fc=t-, 


J  VR^ eikx 


=  —2k  e 


-KS  r+(iK)  r-(iK) 

tR(in)  tR(iK) 


=  -2k  e~KS, 

where  we  have  used  Eqs.(22),  (23)  and  (24).  The  coupling  coefficients  now  take  the  explicit  form 

krl = = Tr  ‘"“4 = ib'"  >•  <52; 


Similarly, 


= wki = it  = it  e'“  )• 


It  is  clear  that  the  coupling  coefficients  consist  of  two  parts:  a  factor  depending  only  upon 
spectral  information  and  waveguide  separation,  and  a  (more  interesting)  contribution  horn  the 
normalization  integral  which  is  dependent  upon  the  actual  geometry  of  the  potential.  To  highlight 
this,  we  will  refer  to  the  normalization  constants  as  shape  factors,  denoted 

F  =  Im  { Res  r_(t«)}.  (54) 

It  is  clear  from  simple  integration  of  Eqs.(42)  and  (43)  that  significant  amounts  of  power  can 
only  be  exchanged  under  conditions  of  phase  matching, 


0L  =  pR  =  $hV  (55) 

For  the  purposes  of  the  optical  interconnect  circuit,  assume  that  the  right-hand  waveguide,  with 
amplitude  A(z)  represents  the  single-mode  branch.  Assuming  the  branch  to  begin  at  some  distance 
z  =  zq  along  the  trunk,  so  that  A{zq)  =  0,  the  coupled-mode  equations  have  the  solutions[14) 

Ai*)  =  B sin  2  ^ 

V  klr  (56) 

B(z)  =  B(zq)  cos  A/?c  z, 
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where 


A/3C  =  \J hrl  *lr  •  (57) 

From  Eqs.(56)  and  (40)  it  is  clear  that  the  field  of  the  composite  trunk/branch  waveguide  closely 
approximates  a  two-mode  system  with  propagation  constants 

/?<+>=  /?  +  A/Jc 

(58) 

/?(-)=  /3-A&. 

Under  conditions  of  complete  power  transfer,  i.e., 

=  *lr,  (59) 

complete  power  exchange  occurs  at  intervals  of 

zl  -  m  =  1,3,5,...  (60) 

2  URL 

measured  from  Zq  along  the  branch  waveguide.  Maximum  power  transfer  is  generally  ensured 
by  employing  identical  waveguides  (not  a  viable  option  in  our  application),  but  it  is  dear  from 
Eqs.(53)  and  (52)  that  equal  coupling  coefficients  can  be  ensured  by  suitable  manipulation  of 
the  normalization  constants  of  the  various  trunk  waveguide  modes  and  the  corresponding  branch 
modes.  For  branches  placed  to  the  right  of  the  trunk,  this  amounts  to 

fjtrunk)  =  ^{branch) .  m  =  1? N .  (61) 

Consequently,  we  will  select  a  set  of  single  mode  branches,  calculate  (gn(branch)  and  suitably 
tailor  the  trunk  waveguide,  using  Eq.(28)  and  selecting  appropriate  values  of  pm  based  upon 
Eq.(61).  In  the  next  section,  we  consider  various  sets  of  branches  and  carry  out  the  trunk  design 
in  accordance  with  these  concepts. 
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6.  DESIGN  EXAMPLES 


The  design  process  consists  of  selecting  suitable  single-mode  branches,  calculating  their 
shape  factors,  and  designing  the  trunk  refractive  index  profile  for  maximum  power  transfer. 
Denoting  the  branch  shape  factors  by  F$n\  this  design  process  imposes  a  condition  on  the  trunk 
through  pm: 


2  Km  Pm 

/i6)  ’ 


(62) 


which  follows  directly  from  Eq.(34)  and  Eq.(61).  It  is  clear  that  the  pm  act  as  the  critical  parameter 
in  matching  waveguides  for  maximum  power  transfer.  We  emphasize  dm  Eq.(62)  is  a  condition 
freely  imposed  upon  the  pm  based  on  the  mode  spectrum  and  the  form  of  the  branch  waveguides. 

At  this  point,  the  goal  of  the  analysis  is  to  determine  to  what  extent  it  is  possible  to  provide 
maximum  coupling  between  the  branch  and  the  trunk  while  adjusting  the  design  parameters  in 
such  a  way  that  die  refractive  index  profiles  are  reasonably  well-behaved. 

A  smooth,  symmetric  trunk  refractive  index  profile  is  dearly  preferable  to  one  with  random 
variations  and  large  gradients.  Using  the  same  mode  spectrum  employed  to  generate  Figs.  3-5, 
we  have  determined  that  a  step-index  branch  design  is  sufficiently  flexible  to  achieve  attractive 
profiles  for  the  trunk  waveguide,  while  maintaining  conditions  of  maximum  power  coupling.  This 
is  an  encouraging  result,  as  step-index  waveguides  are  easy  to  fabricate.  In  all  cases  considered, 
«2  =  3.0,  at  a  wavelength  of  0.9  .urn. 

The  step-index  waveguide  has  been  analyzed  using  standard  methods  (see  [13]).  The  purpose 
of  the  following  analysis  is  to  put  the  step-index  waveguide  into  the  context  of  inverse  scattering 
theory  and  to  demonstrate  its  usefulness  in  the  proposed  interconnect  Consider  a  square  well 
potential  of  width  D  =  2d: 

v(x) 

=  kl  (n2  -  n\)  (- d  <  x  <  d)  (63) 

=  0  elsewhere , 

whose  corresponding  refractive  index  profile  is  a  step-index  planar  waveguide  with  constant  core 
refractive  index  nx.  Defining  the  parameter 

=  \A'2  ~  *u{n2  -  ni }  =  \Aoni  “  (64) 
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wc  can  write  the  Jost  solutions  and  their  derivatives  for  square  well  potential  as 

f-(k,x)  =  aa(k)sinKx  +  ac(k)  cosKx 

(6t 

f'-{k,x)  =  Kaa(k)  cosKx  —  Kac(k)sinKx, 
and 

f+(k,  x)  =  bg(k)  sin  Kx  +  bc(k)  cosKx 

(« 

f+(k,x)  =  Kbt(k)  cosKx  —  Kbc(k)sinKx. 

In  addition,  for  x  <  -d, 

(61 

4(M)  =  -**«-*», 

and  for  x  >  d, 

U(k,x)  =  eikx 

(61 

4(fc,x)=  tfc  e,fcx, 

Continuity  of  the  Jost  solutions  and  their  derivatives  at  these  boundaries  gives  the  coefficients: 


aa(k)  = 


-eifa<{ KsinKd  +  ikcosKd} 


elkd{I(cosKd  —  ikcosKd} 

Oc(fc)  =  - 1 ^ - L, 

b,(k)  =  ~aa(k), 
bc(k)  =  ac(k). 

From  the  first  of  these,  it  is  clear  that  the  eigenvalue  equation  for  the  even  modes  is 


tanKd  = 


The  reflection  coefficient[15]. 


’  W[f_(k,x)J+(k,x)}  ’ 


follows  in  a  straightforward  way.  Since, 


W[f.(k,x),U(k,x)} 

=  K  (bt(k)ac(k)  -  at(k)bc(k)) 

=  2 K  b,(k)ac(k), 
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and 


(73) 


W[f-{-kyx),U{k,x)]  =  2K  (Mc(-fc)  - 


it  follows  that 


r-(k)  = 


_  a,(-k)bc(k)  -  ba(k)ac(-k) 


2  bt(k)ac(k) 


cleaiiy  illustrating  how  the  pole  locations  are  the  respective  even  and  odd  mode  eigenvalue 
equations.  We  are  interested  in  the  fundamental  mode  with  eigenvalue  k\  =  *«i,  for  which  the 


residue  is  given  by 


^  <-<*■»  ■  -r-*u  ™ 

This  expression  can  be  simplified  so  that  the  expression  for  the  coupling  coefficient  takes  the 


K2k2 

l*RL|  =  /?(l  +  «i<0*g(n?-«3) 


e~K,Me2*ld, 


in  exact  agreement  with  the  result  obtained  using  the  standard  method[16]. 

Consider  a  set  of  branches  consisting  of  five  square  wells  of  width  {Dm\m  =  1, 2,  and 
constant  core  refractive  index  =  1, 2, m|.  With  the  eigenvalue  spectrum  preselected, 

the  design  procedure  amounts  to  a  selection  of  branch  core  widths  and  core  refractive  indices 
which  allow  single-mode  operation.  Since  (3^  =  k^nl  +  n^,  and  nffcjj  -/?£,>  0,  the  minimum 


core  refractive  index  is 


■  (^) 


Requiring  the  core  width  to  be  at  least  one  wavelength  nominally  gives 

min{Um}  =  1  fim. 


Table  3  lists  the  results  of  three  sets  of  design  data,  beginning  with  a  set  of  branches  each 
1  fxm  in  width,  their  corresponding  core  retractive  indices,  chosen  so  as  to  satisfy  die  eigenvalue 
equation  for  this  width,  and,  in  the  last  column,  the  corresponding  values  of  pm.  Figure  7  shows 
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the  resulting  refractive  index  profile.  In  the  second  set,  a  similar  pattern  was  followed,  but  for 
larger  cote  widths,  and  the  resulting  trunk  refractive  index  profile  is  plotted  in  Fig.8.  Comparison 
with  the  previous  result  shows  a  greater  shift  of  the  profile  in  the  positive  x  direction,  due  mainly 
to  the  three  pm  which  are  less  than  unity,  in  contrast  to  the  first  case,  where  all  pm  >  1.  Both 
trunks  exhibit  rather  large  index  gradients,  but  are  otherwise  well  behaved. 

The  third  case  is  the  most  interesting.  It  is  clear  from  earlier  discussions  that  if  we  choose 

Pm  at  1  C79) 

for  all  branch  waveguides,  a  smooth,  symmetric  trunk  refractive  index  profile  will  result  Since 
laser  diodes  emit  even  and  odd  field  configurations,  a  symmetric  trunk  refractive  index  profile, 
allowing  for  even  and  odd  guided  modes,  will  result  in  more  efficient  coupling  between  the 
source  and  the  trunk  waveguide.  In  a  somewhat  tedious  but  effective  analysis,  whose  objective 
was  to  satisfy  Eq.(79)  for  all  m,  we  began  by  plotting  a  given  pm  (Eq.(62))  versus  and  Dm. 
Empirically,  it  was  found  that  the  pair  that  satisfied  the  eigenvalue  equation  and  the  condition 
pm  —  1  lay  in  the  vicinity  of  min{n j},  enabling  one  to  narrow  down  the  range  of  Dm  values.  A 
trial  value  of  Dm  was  then  selected,  the  corresponding  found  from  the  eigenvalue  equation. 
Using  this  pair  pm  was  then  checked  for  its  proximity  to  unity.  We  were  satisfied 

to  come  within  3%  of  pm  =  1.  If  required,  the  procedure  can  be  repeated  until  pm  is  sufficiently 
close  to  any  desired  value. 

It  bears  repeating  that  setting  pm  =  1  for  all  modes  merely  guarantees  a  symmetric  trunk 
refractive  index  profile.  The  smoothness  of  the  profile  will  also  depend  upon  the  spectrum  of 
propagation  constants,  as  we  have  seen  in  Figs.3-5.  In  fact,  the  smoothness  of  the  refractive  index 
profile  shown  in  Fig.3  is  a  direct  consequence  of  the  relatively  equal  spacing  of  the  propagation 
constants.  The  more  general  question  of  creating  a  single,  smooth  guiding  region  for  an  arbitrary 
set  of  eigenvalues  and  normalization  constants  is  considered  in  Ref.[7J.  It  is  evident  from  our 
results,  however,  that  the  parameters  governing  the  step-index  branch  waveguides  are  sufficiently 
flexible  to  couple  to  a  large  number  of  possible  trunk  waveguides  designed  using  Darboux 
transformations. 
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7.  CONCLUSIONS 


Guided  wave  optical  interconnects  consisting  of  graded  index  optical  waveguides  were 
designed.  It  was  possible  to  design  an  interconnect  consisting  of  a  multimode  trunk  waveguide 
coupled  to  several  single  mode  branch  waveguides,  each  of  which  delivers  a  selected  mode  to  a 
detector,  and  by  exploiting  the  group  velocity  dispersion  inherent  in  multimode  waveguides,  it 
was  possible  to  select  a  set  of  propagation  constants  such  that  each  of  the  modes  can  be  delivered 
to  its  assigned  detector  simultaneously,  eliminating  clock  skew. 

The  synthesis  of  waveguides  with  prescribed  propagation  constants  is  the  key  to  the  design  of 
this  interconnect.  Consequently,  an  inverse  scattering  algorithm  was  required  to  reconstruct  the 
refractive  index  profile  which  would  support  guided  modes  with  this  preselected  spectrum.  It  was 
determined  that  the  method  of  transformations  provided  a  flexible,  efficient  means  of  generating 
the  multimode  trunk  refractive  index  profiles  suited  to  our  use.  These  profiles  are  continuous  and 
decay  rapidly  in  the  transverse  direction,  making  them  well-suited  to  practical  systems. 

By  manipulating  the  normalization  constants,  it  was  possible  to  take  full  advantage  of  the 
possibilities  of  the  transformation  method.  In  particular,  it  was  possible  to  efficiently  couple  the 
trunk  waveguide  and  each  of  the  branch  waveguides,  despite  the  fact  that  the  trunk  and  branches 
consisted,  in  general,  of  different  refractive  index  profiles.  This  analysis  resulted  in  a  formulation 
of  waveguide  coupling  coefficients  in  terms  of  the  scattering  data  pertaining  to  the  corresponding 
potentials.  It  is  emphasized  that  this  formulation  is  completely  general  and  applicable  to  any 
waveguide  systems  in  which  the  weak  coupling  approximation  is  valid. 

In  addition,  it  was  found  that  proper  manipulation  of  the  normalization  constants  guaranteed 
trunk  refractive  index  profiles  which  were  symmetric  and,  under  certain  circumstances,  free  from 
large  index  gradients. 

Directions  for  future  work  include  analyzing  the  sensitivity  of  the  refractive  index  profiles 
to  variations  in  the  propagation  constants,  and  an  in-depth  analysis  of  allowed  variations  in  chip 
placement  within  the  prescribed  wafer  area. 
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pm(xkon2),pm. 

Figure  3 

Figure  4 

Figure  5 

0s,  Ps 

1.09752,  1 

1.05169,  1 

1.09752,  547 

04,  Pa 

1.07611,  l 

1.05269,  1 

1.07611,  0.003 

(h,P3 

1.05369,  1 

1.05369,  1 

1.05369,  0.01 

02, Pi 

1.03202,  1 

1.05469.  1 

1.03202,  96 

Pi 

1.01034,  1 

1.05459,  1 

1.01034,  235 

Table  1.  Data  characterizing  the  refractive  index  profiles  of  Figs.3-5. 
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SMOOTH  PROFILE,  (Figure  3),  Direct  method  comparison 

An  (x*on2) 

(3m  (x  A^n2),  direct  method 

1.09752 

1.09752 

1.07611 

1.07614 

1.05369 

1.05377 

1.03202 

1.03218 

1.01034 

1.01054 

Table  2.  Propagation  constants  for  smooth  profile.  Second  column 
gives  results  of  finite  difference  method  applied  to  potential  in  Fig.3. 


STEP  INDEX  BRANCH  DATA 


An  (xfc0n2) 

d  (=  D/2),  fim 

Til 

Pm 

Fig.7 

1.09752 

0.5 

3.31 

30.77 

1.07611 

0.5 

3.25 

99.81 

1.05369 

0.5 

3.18 

107.95 

1.3202 

0.5 

3.11 

7.83 

1.01034 

0.5 

3.04 

8.83 

|  Fig.8 

1.09752 

1.0 

3.30 

5.89 

1.07611 

1.0 

323 

226 

1.05369 

1.0 

3.17 

0.54 

1.3202 

1.0 

3.10 

0.01 

1.01034 

1.0 

3.03 

0.003 

|  Fig.3,  approx.  | 

1.09752 

0.648 

3.30613 

1.00901 

1.07611 

0.664 

3.24113 

1.00456 

1.05369 

0.947 

3.16782 

0.984028 

1.03202 

1.1 

3.10102 

0.97903 

1.01034 

1.4 

3.03384 

1.00375 

Tabic  3.  Data  for  step-index  branches  corresponding  to  trunks  in  Figs.7,  8,  and  3. 
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Figure  Captions 


1.  The  optical  interconnect  Refractive  index  profiles  are  designed  so  that  a  pulse  launched 
from  point  S  teaches  each  of  the  points  P(i)-..P(n)  simultaneously. 

2.  One  dimensional  planar  waveguide  with  variable  core  refractive  index  n(x)  surrounded  by 
cladding  layers  of  constant  refractive  index  n^. 

3.  The  smooth,  symmetric  trunk  refractive  index  profile  resulting  from  evenly  spaced  eigen¬ 
values  and  pm  =  1  for  all  m. 

4.  Think  refractive  index  profile  resulting  from  five-fold  near  degeneracy.  Symmetry  is  retained 
since  pm  =  1  for  all  m. 

5.  Trunk  refractive  index  profile  with  same  spectrum  as  in  Fig. 3.  The  pm  are  varied  as  indicated 
in  Table  1. 

6.  Weakly-coupled  waveguides  used  to  model  the  coupling  interaction. 

7.  Think  refractive  index  profile  for  step  index  branch  waveguides  of  width  1  pm. 

8.  Trunk  refractive  index  profile  for  step  index  branch  waveguides  of  width  2  pm. 
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igure  7.  Mills  f.  Tamil 
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Abstract 


Within  the  context  of  weak  coupling  theory,  we  derive  representations  of  the  coupling 
coefficients  between  neighboring  waveguides  by  representing  the  field-dependent  inter¬ 
action  integrals  by  algebraic  expressions  involving  scattering  data  and  we  illustrate  the 
contexts  in  which  scattering  theory  can  make  a  viable  alternative  to  existing  formulation 
of  the  waveguide  coupling  problem. 
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I.  Introduction 


Coupling  between  waveguides  in  a  multilayer  structure  is  the  cornerstone  of  optical 
spatial  switching.  This  form  of  coupling,  which  arises  when  the  evanescent  fields  of 
one  waveguide  perturb  its  neighbor,  can  be  analyzed  by  several  methods.  Traditionally, 
the  most  popular  approaches  have  been  a  weak  coupling  perturbation  analysis[l]  or  the 
analysis  of  local  normal  modes[2].  The  problem  continues  to  generate  considerable  in¬ 
terest,  as  evidenced  by  recent  work  formulating  variational  methods  and  finite-difference 
schemes[3].  The  results  presented  in  this  paper  were  motivated  by  our  recent  interest 
in  analysis  of  optical  waveguides  using  scattering  data  (i.e.,  eigenvalues  of  the  bound 
modes,  reflection  and  transmission  coefficients). 

Specifically,  this  paper  shows  that  the  traditional  weak-coupling  analysis  of  interact¬ 
ing  waveguides  can  be  reformulated  in  the  language  of  scattering  theory.  We  show  that 
the  coupling  coefficients  describing  the  interaction  of  two  neighboring  waveguides  have 
straightforward  representations  in  terms  of  their  scattering  data,  eliminating  the  need  to 
explicitly  calculate  the  field-dependent  interaction  integrals  by  representing  these  inte¬ 
grals  with  straightforward  algebraic  expressions  involving  the  guided-mode  propagation 
constant  and  the  residue  of  the  reflection  coefficient  In  this  paper  no  attempt  is  made  to 
reformulate  the  mathematics  of  scattering  theory,  but  rather  to  identify  existing  aspects 
this  theory  which  are  useful  when  applying  transverse  coupling  to  waveguide  design, 
and  illustrating  the  contexts  in  which  scattering  theory  can  make  a  viable  alternative  to 
existing  methods. 


II.  Waveguide  Model 

Figure  1  shows  a  planar  waveguide  consisting  of  two  coupled  graded-index  (GRIN) 
guiding  regions.  For  the  moment,  consider  a  single  planar  graded-index  waveguide 
consisting  of  an  inhomogeneous  core  with  a  varying  refractive  index  n(x),  surrounded 
by  two  cladding  layers  of  constant  refractive  index  n2.  To  simplify  the  analysis  we 
assume  that  each  guiding  region  is  infinite  in  the  y  direction  and  supports  a  single  y— 
polarized  TE  mode  of  the  form 


£,(*,»,«)  s  £,(*)<*•«-" 1  ,  (1) 

where  z  is  the  direction  of  propagation,  u  is  the  frequency,  /?  the  longitudinal  propagation 
constant.  It  has  been  assumed  that  the  waveguide  is  infinite  in  extent  along  the  y 
axis  Here,  is  the  free  space  wavenumber.  The  field  Ey(x)  is  defined  by  the  scalar 
differential  equation 


+  [ifcoM*) W  =  0. 


(2) 
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This  equation  can  take  the  form  of  a  Schrodinger  equation  which  is  particularly  well 
suited  to  analysis  using  scattering  data.  Defining  the  complex  transverse  propagation 
constant  k  (=  kr  -f  ikt)  as 


k2  =  kln\  -  ft2  (3) 

brings  Eq.(2)  into  the  Schrodinger  form 

^+[*2-„(*)]E,  =  0  (4) 

whose  potential 

v(x)  =  fcjj  [n22  -  n  2(x)]  (5) 

varies  across  the  waveguide  core  and  vanishes  in  the  cladding.  Equation  (5)  clearly 
illustrates  how  the  depth  of  the  potential  may  be  varied  either  by  changing  the  wavelength, 
altering  the  refractive  index  profile,  or  both.  In  this  scheme  the  mode  cutoff  condition, 
ft  =  k0n 2,  is  obtained  when 


k  =  0.  (6) 

llie  discrete  set  of  guided  modes,  characterized  by  fcon2  <  ft  <  k0ni,  or  equivalently  by 

0  <  Im  k  <  Im  yjn\  —  n^j  ,  is  represented  by  points  along  the  positive  imaginary 

axis  of  the  complex  k  plane.  In  scattering  theory,  the  guided  modes  are  termed  bound 
states,  distinguished  by  their  discrete  eigenvalues  k.  As  the  fundamental  mode  of  a 
planar  waveguide  is  TE,  Eq.(4)  is  sufficient  to  describe  the  bound  mode  in  a  single-mode 
waveguide. 

A.  Scattering  Coefficients  and  Jost  Solutions  Scattering  theory  (direct  and  inverse)  is 
concerned  with  the  relationship  between  a  Schrodinger  potential  t>(x)  and  its  associated 
scattering  data  (i.e.,  reflection  and  transmission  coefficients).  A  plane  wave  e+ikx  incident 
on  the  potential  from  x  =  — oo,  will  give  rise  to  a  reflected  portion  taking  the  form, 

r_{k)e~ikx  (7) 

as  x  — >  —  oo,  as  well  as  a  transmitted  wave, 

t-{k)e+lkx  (8) 

as  x  — ►  oo.  An  alternative  viewpoint  is  provided  by  the  coefficients  r+{k )  and  t+(k) 
which  define  relected  and  transmitted  portions  of  a  plane  wave  incident  from  x  =  +oo. 
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The  Schrodinger  equation  admits  a  pair  of  Jost  solutions,  denoted '/+(£,  x)  and 
f-(k,  x),  defined  according  to  their  asymptotic  behavior: 

lim  f+(k,x)e~tkx  =  Iim  f-(k,  x)e+‘*r  =  1.  (9) 

i— >-4-oo  x— »— oo 

The  pairs  {/+(&,  x), /+(-&,  x)}  and  {/_(&,  x), /_(-&,  x)}  comprise  sets  of  linearly 
independent  solutions  to  the  Schrodinger  equation,  allowing  construction  of  the  linear 
combinations^]: 


f±{k,x)  =  ^/T(-fc,x)  +  ^®/T(fc,x),  (10) 

The  Wronskian,  defined  as  W[f,g]  =  f  g'  —  g  /',  (the  prime  denoting  differentiation  with 
respect  to  the  coordinate),  provides  a  set  of  relations, 

TM-T^)  =  w^h^k^  <■" 

so  that  f_(fc)  =  t+(k)  =  t(k),  a  result  which  is  a  direct  consequence  of  the  asymptotic 
behavior  stipulated  in  Eq.(l-2).  In  addition, 

“if  =  02) 

follow  from  Eqs.(10)  and  (11). 

During  the  course  of  this  analysis,  it  is  useful  to  shift  potentials  along  the  axis.  The 
scattering  data  changes  in  a  controlled  way  under  a  shift  Consider  a  potential  v(x)  with 
Jost  solutions  denoted  f±(k,x)  and  corresponding  scattering  data  r±(k),t(k).  It  is  clear 
that  the  shifted  potential  v(x  —  d)  has  a  Jost  solution  of  the  form 

f+{k,x)  =  elkdf+(k,x  -  d).  (13) 

Using  an  overbar  to  denote  the  scattering  coefficients  of  the  shifted  potential,  it  can  be 
shown  that  the  reflection  coefficients  associated  with  the  translated  potential  are  related 
to  the  original  data  by  a  simple  phase  shift 

r.(k)  =  e+2M  r,(k), 
f+(k)  =  e-2Mr+(k), 

while  the  transmission  coefficient  is  unaltered: 

t{k)  =  t{k).  (15) 
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B.  Guided  Modes 

At  values  of  k  (=  ia)  such  that 


(16) 


that  is,  the  bound  state  eigenvalues  correspond  to  the  poles  of  t(k)  which  lie  on  the 
positive  Im  k  axis.  The  Jost  solutions  exhibit  the  asymptotic  behavior 


Ebyound{x )  ~  e*ai,  x  -►  ±oo  (a  >  0), 

This  implies 

...  r_(za)  ,  , .  \ 

U(ia,x)  =  f-(ia,x), 

f-{ia,x)  =  7 f+(ia,x), 

resulting  in  the  following  useful  relation: 

r-fo)  M»a)  j 

t-(ia)  f+(ta) 

The  corresponding  normalized  guided  mode  fields  are  then 

Ey0Und(x)  =  c+  U(ia,  x)  =  c_  f-(ia,  x), 


(17) 


(18) 


(19) 


(20) 


where  the  c±  are  arbitrary  constants.  There  is  a  one-to-one  correspondence  between  the 
bound  states  of  the  quantum  mechanics  picture  and  these  guided  modes,  and  we  will  use 
these  terms  interchangeably. 


m.  Transverse  Coupling 

In  this  section  we  review  the  salient  features  of  waveguide  interactions  in  the  weak 
coupling  approximation  and  reformulate  the  problem  in  terms  of  scattering  data. 

Figure  2  shows  two  neighboring  (non -overlapping)  potentials,  separated  by  a  distance 

s: 

s  =  dR-  di.  (dL  <  0)  (21) 

Each  waveguide  is  assumed  to  have  a  graded-index  core  with  refractive  index  profiles 
ni(x)  and  tir(x): 

n\{x)  =  nf  +  A  n2L(x) 
n\{x)  =  nf  +  An^(x). 
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We  will  assume  that  each  waveguide  separately  supports  y—  polarized  TE  modal  fields 
El(x)  and  ER(x)  with  propagation  constants  Pi  and  PR,  respectively.  The  interaction 
between  the  two  fields  will  be  represented  by  a  z—  dependent  linear  combination  of  the 
individual  waveguide  modes: 

S(x ,  2,  t)  =  A(z)  ER{x)e^wi-^z)>  +  B{z)  EL{x)e*ut-fiL*\  (23) 


the  exact  form  of  the  z—  dependent  weighting  coefficients  A(z)  and  B(z)  being  a  function 
of  the  interaction  strength  induced  by  the  transverse  coupling.  The  governing  equation 
for  the  field  in  Eq.(23), 


d2s 

dx 2  ^ 


d2S  w2  2 

dz2  n  ~0, 


where  n2(x)  is  the  refractive  index  of  the  composite  structure. 


n2(x)  =  n2  4-  An\(x)  -f  A n2R(x), 


(24) 

(25) 


dictates  the  coupling  analysis  which  carries  with  it  three  explicit  assumptions: 

i)  Each  waveguide  individually  supports  single  mode  with  propagation  constants  Pi 
and  pR. 

ii)  The  coupling  is  weak,  i.e.. 


<PA 

dA 

<PB 

„  dB 

dz 2 

« 

0Ri; 

dz 2 

<< 

(26) 


and  iii), 

+oo  +oo 

J  El(x)  Ei(x)  dx  «  J  E2l  r(x)  dx. 

— oo  —  oo 


(27) 


Ignoring  the  second  derivatives  of  A{z)  and  B(z)  the  wave  equation  reduces  to  a  set 
of  coupled  first-order  differential  equation  for  the  z-dependent  coefficients: 


—  =  i  krl  B(z)e  X^L  0r)z  +i  krrA(z), 

J  D 

—  =  i  kiR  A(z)e'^L~0R)z  +  i  Km  B(z). 
where  the  coupling  coefficients 


(28) 


URL 


KRR 


lRL 

2  PrNrr 
Irr 

^PrNrr 


(29) 
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arc  functions  of  the  interaction  and  normalization  integrals: 


+OO 


IRL  =  J  Er(x)vr(x)  EL(x)dx 

— oo 
+oo 

Irr  =  J  Er(x)  vr(x)  Er(x)  dx 

— oo 
+oo 

Nrr  =  J  Er(x)  dx. 


(30) 


The  coefficients  klr  and  kh  follow  by  interchanging  L  and  R.  The  self  coupling 
coefficients  kh  and  krr  represent  small  corrections  to  the  propagation  constants  and 
are  usually  ignored. 

In  the  phase  matched  case  {0i  =  0r  =  0),  the  solutions  to  Eq.(28)  are 


A(z)  =  ^  A(zo)  cos/S.0  z  +  B(zo)  smA/3z| 

B(z)  =  |b(z0)  cosA0  z  +  iyl- A(zo)  smA/?z| 


(31) 


where 


A/3  =  y/*RL  *LR  ■ 


(32) 


Given  the  initial  condition  A(zq)  =  0,  the  solutions  become 

A(z)  =  i  B(zo)  sin  A/?  z 
B(z)  =  B(zo)  cosA/3  z, 

provided 


(33) 


URL  =  KLR  =  K,  (34) 

so  that  complete  power  transfer  occurs  at  intervals  of  (tt/2)  A 0  along  the  coupled  length 
of  the  waveguides.  Substituting  these  expressions  for  A(z)  and  B(z)  into  Eq.(23) 
indicates  that  the  total  electric  field  consists  of  an  approximate  linear  superposition  of 
modes  with  propagation  constants. 


0+  =  0  +  A0 

0~  =  0  -  A/?. 


(35) 


The  interaction  and  normalization  integrals  are  the  principal  calculational  hurdle 
associated  with  the  weak  coupling  model.  The  latter,  as  previously  shown,  have  a 
convenient  representa'.ion  in  terms  of  the  scattering  data.  In  the  next  paragraphs  we 
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outline  established  scattering  relationships  which  are  useful  in  representing  the  interaction 
integrals. 

Consider  two  different  potentials  v(x)  and  v(x)  (These  are  not,  in  general,  the 
potentials  describing  the  two  individual  waveguides)  [4].  In  the  limit  as  x  — >  —  oo, 
it  follows  that 

W  [/+(*,  x),  /+(*,  x)]  =  -J±-{r_(k)  -  f(i)],  (36) 

while  in  the  limit  as  x  — >  oo,  this  Wronskian  vanishes  since  f+(k,x)  =  f+(k,x)  ~  e'hz. 
Now  consider  the  derivative, 

-^W  [/+(*:,  x),  /+(fc,  x)j  =  [u(x)  -  u(x)]  f+(k ,  x)/+(fc,  x).  (37) 

Integrating  this  expression  yields 

+oo 

/  ^w[Mk<*)J+(k,x)]dx  =  -w[Mk,x),Mk,x)]  (38) 

— OO 

providing  a  convenient  integral  relation  for  the  Wronskian: 

+oo 

/2  ik 

[u(x)  -0(x)]  f+(k,x)f+{k,x)dx  =  [r_(fc)  -f_(fc)].  (39) 

—  OO 


Using  similar  steps,  a  companion  expression  can  be  derived: 


+oo 

J  [u(x)  —  v(x)]  /_(&,  x)f-(k,  x)  dx  = 

—  OO 


2  ik 

t(k)t(k) 


[r+{k)  -  f+(k)]. 


(40) 


If  v(x)  =  0,  Eqs.(39)  and  (40)  reduce  to  the  useful  form, 

+oo 

J  f±{k,x)v{x)e±,kx  dx  =  ^  r^(k).  (41) 

—  OO 


Writing  the  guided  mode  electric  fields  of  the  right-  and  left-hand  waveguides  in 
terms  of  the  Jost  solutions  for  the  respective  potentials  v°R(x)  and  v°L(x)  gives 


ER(x)  =  f°_R(ia,x  -  dR) 
EL(x)  =  f°+L(ia,x-dL), 


(42) 
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where  a  bound  state  with  eigenvalue  k  =  ia  (a  >  0)  is  assumed.  With  the  help  of 
Eqs.(19)  and  (41)  we  find: 

+oo 

Ilr  =  e~adR  j  flL(k,x-  dL)  vL(x)e-ikxdx\k=ta  =  -2a  e"",  (43) 

—  OO 


The  normalization  integral  Nil  is  simply. 


4-oo 

Nll  =  J  [/JL(ta,x  -  dz,)]  dx 

— OO 


i 

Res  r+L(ia)  ’ 


(44) 


Defining  the  shape  factors 


F3  =  Im  j/Ees  r!^(ia)  j 
F+  =  Im  |  Res  (ia)  | , 


gives 


and  by  the  same  token. 


KLR  = 


krl  = 


Ilr 


20  Nll  0 
Irl 


~a  e~as  Fi, 


20  Nrr  0 


^  e~as  Ff . 


(45) 

(46) 

(47) 


These  results  indicate  that  when  placed  in  the  context  of  scattering  theory,  the  directional 
coupling  coefficients  take  on  a  particularly  simple  form  consisting  of  a  portion  which 
is  dependent  solely  upon  the  eigenvalue  and  the  waveguide  separation,  multiplied  by  a 
factor  whose  value  is  dependent  upon  the  inherent  shape  of  the  potential.  The  condition 
for  complete  power  transfer  takes  the  particularly  simple  form. 


(48) 


One  form  of  this  condition,  which  is  likely  to  be  encountered  in  practice,  is  simply 

r«(k)  =  rj(*),  (49) 

which  implies 


vl(x)  =  vr(-x ), 


(50) 


and  the  intuitively  appealing  conclusion  that  a  waveguide  is  coupled  with  100%  efficiency 
to  its  “mirror  image”. 
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IV.  Design  Examples 


A.  Step-Index  Waveguides 

Coupling  in  step-index  waveguides  has  been  analyzed  using  standard  methods  (see 
[5]),  therefore,  it  can  serve  as  a  check  for  the  coupled  mode  formalism  we  have  derived. 
Consider  a  square  well  potential  of  width  D  =  2 do: 


v(x)  = 


&o  (n2  _  ni)  (_^o  <  x  <  do) 
0  elsewhere , 


whose  corresponding  refractive  index  profile  is  a  step-index  planar  waveguide  with 
constant  core  refractive  index  ni.  Defining  the  parameter 

K  E  y/p-i =  yjq n\  -  0\  (52) 

we  can  write  the  Jost  solutions  for  square  well  potential  as 

a3(k)  sin  Kx  +  ac{k)  cosKx  (—do<x<do) 

/-(fc’l)=  «-»*  (,<-*),  <s3> 


b3(k)  sin  Kx  +  bc{k)  cosKx  (—do<x<do) 
f+{k,x)=  e+ikx  (x>dQy 


Continuity  of  the  Jost  solutions  and  their  derivatives  at  these  boundaries  gives  the 
coefficients: 

— etkd°{KsinKdo  +  ikcosKdo} 
a3{k)  = - — - , 


e'kd°{KcosKdo  —  iksinKdo} 

ac(k)  = - - - ,  y 

b3(k)  =  —a3(k), 
bc(k)  =  ac(k). 

The  reflection  coefficient  (from  Eqs.(l  1  and  (12)),follows  in  a  straightforward  way: 

_  as{-k)bc(k)  -  b3(k)ac{—k) 
r'[k)  ~  2 b3{k)ac(k) 

whose  pole  locations  lead  to  the  familiar  eigenvalue  equations 


tan  I\d(i  = 


-IK 

—  (even) 

r  r 

4  {odd) 
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(57) 


This  result  is  applicable  to  both  single-  and  multimode  waveguides.  We  are  interested  in 
the  fundamental  mode  with  eigenvalue  k\  =  ia,  for  which  the  residue  is  given  by 


_  _!  „-2ikd0  _ 

2  d 


Res  {r_(ia)}  =  —re 


Ktan  I<do  —  ik 


tanKdo  -f  ik} 

With  the  help  of  the  Eqs.(47)  and  (57),  the  coupling  coefficient  takes  the  form 


(58) 


I«kz,I  = 


AV 


e-ase2ad0 


/?(1  +  ad)  kl(n\  -  nl) 
in  exact  agreement  with  the  result  obtained  using  the  standard  method[5]. 


(59) 


B.  Depressed  Cladding  Waveguides 

The  foregoing  result  puts  scattering  theory  in  direct  contact  with  established  results, 
but  provides  little  motivation  to  apply  scattering  theory  as  opposed  to  the  conventional 
techniques,  due  largely  to  the  fact  that  the  guided-mode  fields  have  straightforward 
representations  and  the  interaction  and  normalization  integrals  can  be  readily  calculated. 
As  refractive  index  profiles  become  more  complicated,  the  need  for  an  alternative  method 
becomes  more  compelling.  Jordan  and  Lakshmanasamy[6]  designed  high  V-number 
single-mode  planar  waveguides  using  a  rational  reflection  coefficient  of  the  form 


r-(k) 


—  k\  k2k3 


{k  -  ki)(k  -  k2)(k  -  k3y 


(60) 


which  yields  a  single  bound  mode  eigenvalue  at  k3  =  ia,  and  two  poles  ki—  —  cj  —  ic2 
and  fc]=  ci  —  ic2  in  the  lower  half  of  the  complex  Iplane  which  represent  tunneling  leaky 
waves.  The  authors  showed  that  the  Gelfand-Levitan  reconstruction  technique  results  in 
a  corresponding  potential 


v(x)  =  2 


^  -  q(x)A-'(x)A'(x) 


A  1(x)7/, 


(61) 


where  a  and  7  are  the  row  vectors 

a  =  [1  1  eVlX  e~VlX  e*x  e"’** 

7  =  [00  0  0  0  -  a(c[+^)  ], 

and  A(x)  is  a  6  x  6  matrix  whose  elements  are  listed  in  Appendix  1.  The  parameters 
are  defined, 


(62) 


m  =  [(<r  +  p)/ 2]1/2 
m  =  {(<r  -  p)/2]1/2 
<T  =  a2  +  2cl~  2c] 

P  ~  [(°2  —  4c^)  (a2  +  4c] )] >/2. 


(63) 
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Some  restrictions  apply  to  the  relative  location  of  the  poles  a,  c\  and  c2  brought  about 
by  requiring  a  real  potential.  Specifically, 

0  <  C2  <  —■  (64) 

places  a  lower  bound  on  c2  and 

<r  >  p  (65) 

etches  an  upper  limit  on  c\  and  c2.  This  condition  is  identical  to  the  one  derived  in  [6] 
based  upon  the  conservation  of  energy  condition, 

|r(/b)|2  <  1  all  Re  k.  (66) 

Each  of  the  three  refractive  index  profiles  illustrated  in  Fig.3  propagates  a  single  mode 
with  propagation  constant 


/?=  1.01034  k0n2. 


(67) 


The  depressed  portion  of  the  refractive  index,  characterized  by  a  portion  of  the  profile 
dipping  below  the  nominal  AlGaAs  cladding  value  n2  =  3.0,  is  most  clearly  evidenced 
as  the  poles  for  the  tunneling  leaky  modes  are  moved  farther  from  the  lower  Im  k  axis. 

The  residue  at  the  pole  representing  the  bound  mode  is  easily  found  to  be 


Res  r_(ia)  =  ia 


(4  +  4) 

(4  •+  (a  +  c2)2) 


a2  +  72 

a2  +  (1  +  7)2 


(68) 


where  ci  =  a  a  and  c2  =7  a.  In  Fig.  4,  the  shape  factor  is  plotted  as  a  function  of 
ci,  showing  a  monotonic  increase  (for  a  given  c2)  as  the  poles  are  moved  out  into  the 
complex  plane.  Based  on  the  form  of  the  refractive  index  profiles  themselves,  this  result 
is  expected,  as  a  decrease  in  ci  is  associated  with  translation  of  the  optical  channel  along 
the  positive  x  axis. 

Figure  3  suggests  that  for  small  values  of  ci  and  c2,  the  refractive  index  profile 
approximates  a  seek2  x  form,  suitably  scaled  and  translated  a  finite  distance  along  the 
positive  1  axis.  This  is  indeed  the  case,  and  such  profiles  (developed  from  a  slightly 
different  perspective)  are  taken  up  in  the  next  subsection. 
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C.  Truncated  Refractive  Index  Profile  Waveguides 

We  now  consider  the  family  of  single-mode  refractive  index  profiles  based  on 
truncated  versions  of  the  potential 

v3(x)  =  —  4i2  sech2  by/2  x,  (69) 

parameterized  by  the  positive  scaling  constant  b.  Potentials  of  this  form  are  single-mode 
with  a  bound  state  eigenvalue  k  =  i  by/2.  Equation  (69)  is  representative  of  a  smooth 
function  which  decays  relatively  rapidly  for  large  |xj  making  it  a  suitable  refractive 
index  profile  for  waveguide  design.  For  the  purposes  of  this  paper,  a  truncation  is  a 
discontinuity  imposed  upon  a  smooth  potential,  at  a  point  x  —  x\  such  that 

v(x)  =  0  (x  <  xi).  (70) 

Clearly,  this  creates  a  cladding  region  of  constant  refractive  index 

n(x)  =  n.2  (x  <  xi).  (71) 


In  previous  work  we  completed  an  extensive  analysis  of  these  potentials  from  the 
standpoint  of  scattering  theory,  including  the  effects  of  truncations  of  the  potentials  to 
model  core-cover  interfaces,  considering  both  single-  and  multimode  waveguides[7].  In 
the  present  paper,  we  extend  our  analysis  to  include  the  effects  of  truncations  upon  the 
coupling  coefficient  Although  we  restrict  ourselves  to  single-mode  waveguides  arising 
from  potentials  of  the  form  in  Eq.(69),  coupling  between  modes  in  multimode  waveguides 
follows  a  similar  analysis. 

We  have  previously  shown  that  the  transmission  coefficient  may  be  written  in  terms 
of  the  Jost  solution  /+(/:,  x)  of  the  corresponding  untruncated  potential.  In  the  case  of  a 
single  truncation  at  the  point  x  =  x\,  the  transmission  coefficient  takes  the  particularly 
simple  form: 


tT(k)  =  2ik  eikx'  [f\.(k,x1)  +  ikf+(k,x1)]  \  (72) 


whose  poles  provide  the  eigenvalues  (and  corresponding  propagation  constants)  as  a 
function  of  x\.  Here  the  prime  denotes  differentiation  with  respect  to  the  spatial 
coordinate.  Up  to  this  point  the  results  are  completely  general. 

The  Jost  solutions  for  potentials  v.,(x)  of  Eq.(69)  take  the  form. 


f+{k,x) 


ikz 


ik  —  6\/2  tank  b\/2  x 
ik  —  by/2. 


(73) 
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Used  in  conjunction  with  Eq.(72),  the  Jost  solution  provides  an  analytic  expression  for 
the  bound  state  eigenvalue  k  , 


k  =  i  — 
2 


—b\/2  tank  b\/ 2  xj  +  b\pl  \J \  +  seek 2  6\/2 


Xl 


(1A\ 


The  ratio 


r^(fc)  e‘*Zl  ,  .  , 

W{k)  =  aT  t^/+(^Xi)  -  /+(fc^i)] 

combined  with  Eqs.(72)  and  (73)  gives  the  reflection  coefficient 

n  _  -be2ikx'  seek2  by/2  xx 
r~U  (A:  ~  kp)  (k  -  kn)  ’ 


(75) 


(76) 


where  A:p  and  kn  lie  on  the  positive  and  negative  imaginary  axes,  respectively,  taking 
the  values: 


l-bV2  \  ~tl±J  (l  +sech2bV2xl)  ^  tl  =  tanhby/ixx.  (77) 


When  Eq.(72)  is  reduced,  the  transmission  coefficient  has  the  simple  form. 


tT(k)  = 


2k  (k  +  i  by/2) 

( k  kp)(k  kn ) 


(78) 


The  shape  factor  of  the  branch  can  be  conveniently  expressed  as  a  function  of  the 
truncation  point, 


pT  _  be2,kpXl  seek2  b\/2x\ 
^2(l  +  seek 2  b\/2x\ ) 


(79) 


In  Fig.  5  we  plot  the  shape  factor,  along  with  2  kp/i,  as  a  function  of  the  truncation  point 
xj.  For  comparison,  we  have  included  the  magnitude  of  the  area  under  the  potential: 


—46  r 

— -pr  l—tanhbV2xi 
v2 


(80) 


which  is  also  a  monotonically  decreasing  function  of  xx.  It  is  interesting  to  note  that  the 
decrease  in  the  shape  factor  more  closely  parallels  the  behavior  of  the  area  for  a  larger 
interval  of  x\  than  it  does  the  eigenvalue  itself. 

In  section  III  we  emphasized  that  the  well-known  coupled-mode  electric  field  is 
effectively  a  two-mode  solution  for  the  composite  double  -  well  system  representing  the 
coupled  waveguides.  For  two  reasons,  it  is  appropriate  to  follow  up  on  the  implications  of 
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Eqs.  (32).  First,  scattering  theory  provides  a  straightforward  way  to  evaluate  eigenvalues 
of  a  composite  structure  consisting  of  two  non-overlapping  potentials,  and  second,  it 
provides  further  verification  that  scattering  analysis  of  the  coupled  mode  problem  is 
consistent  with  known  results. 

Starting  from  first  principles,  it  is  straightforward  to  show  that  the  transmission 
coefficient  of  the  composite  potential  (see  Fig. 2)  takes  the  form[8], 

oo  m  /  1  \ 

m  =  tL(k)tR(k)  £  (r?(t)  r i(k))  =  tL(k)tR(k)  (1_rJiWrtw)  •  (8» 


The  scattering  data  in  this  expression  applies  to  the  potentials  vi(x)  and  vR(x).  From 
this  point  we  will  reduce  this  general  result  to  encompass  the  special  case  of  two  “mirror 
image”  single-mode  truncated  potentials  separated  by  a  distance  2d  (i.e.,  dR  =  —di  =  d) 
for  which  the  scattering  coefficients  take  the  form  of  fractions  made  up  of  arbitrary 
k  dependent  functions,  the  numerator  and  denominator  denoted  with  the  appropriate 
subscripts  n  and  d: 


rO:(k)  =  r+(k)  =  e 


.2  kd 


rn{k) 

rd{k) ’ 


and 


tLM  =  tRM  =  Wy 


so  that  the  composite  transmission  coefficient  can  be  written, 

«(*) 


m  = 


rj(*)  -  fi(k) 


(82) 

(83) 

(84) 


It  is  clear  that  for  sufficiently  large  separations,  the  transmission  coefficient  will  exhibit 
two  closely-spaced  bound  state  eigenvalues  lying  close  to  the  original  single  eigenvalue. 
As  we  mentioned,  the  corresponding  propagation  constants, 


/?+  =  P  +  A/?+ 

r  =  /?-A/r, 


(85) 


which  are  roots  of  the  denominator  of  Eq.(84),  provide  an  approximation  to  the  coupling 
coefficient(see  [1]) 


k  ~  A/3+  ~  A/3  , 


(86) 


provided  that  the  waveguides  are  weakly  coupled. 

Consider  the  truncated  single-mode  potential  of  Eq.(69), 


—  462  sech2  6\/2  x 

vR(x)  = 


x  >  0 
x  <  0. 


(87) 
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whose  reflection  coefficient 


'-<*>  =  FT6-  (88) 

follows  directly  from  Eq.(76)  with  x\  =  0.  There  is  a  single  bound  state  eigenvalue  at 
k  =  ib.  The  composite  structure,  consisting  of  this  potential  shifted  a  distance  d  along 
the  +x  direction,  and  its  mirror  image,  each  with  reflection  coefficients  (see  Eq.(14) 

r£(*)  =  r?(*)  =  e*“1^,  (89) 

has  eigenvalues  k ±  =  fd±2  —  kfol  given  by  the  roots  of  the  denominator  of  Eq.(84). 

An  approximation  to  the  shape  factor  of  the  single-mode  potential  is  found  by 
inverting  Eq.(47): 

F  ~  ^  eba  Afl+  ~  ^  ebs  A/T,  (90) 

b  b 

where  s  =  2d.  As  the  separation  is  increased,  a  convergence  of  /3+  and  /?-  towards  f3  is 
expected,  leading  to  better  approximations  to  the  shape  factor. 

In  table  1  we  list  the  eigenvalues  and  corresponding  approximate  values  of  the  shape 
factor  (from  Eq.(90))  against  the  exact  value 

Im  Res  r^(ib)  =  ^  =  1.889  x  106.  (91) 

(We  have  taken  b  =  3.778  x  106).  As  d  is  increased,  the  expected  rapid  convergence  to 
the  correct  shape  factor  is  readily  apparent  Physically,  this  is  the  result  of  the  increasing 
accuracy  of  the  weak  coupling  model  as  the  waveguides  are  separated. 

Aside  from  the  general  analytic  interest  of  this  approximation,  it  may  be  advantageous 
to  apply  it  in  situations  where  the  residue  of  the  reflection  coefficient  is  a  complicated 
or  difficult  to  calculate.  In  our  experience,  this  is  often  the  case  for  refractive  index 
profiles  incorporating  two  truncations  (to  simulate  two  cladding  regions),  for  which  the 
residues  undergo  extremely  rapid  variations  in  the  vicinity  of  the  bound  state  eigenvalues. 
Certainly,  further  work  is  needed  in  this  area. 

V.  Conclusions 

Coupling  in  multilayer  waveguide  structures  is  studied  here  using  scattering  tech¬ 
niques.  As  inverse  methods  find  wider  applications  in  waveguide  design,  the  inverse 
scattering  representation  of  transverse  coupled  modes  developed  here  will  be  useful  in 
the  design  of  multilayer  devices.  By  replacing  the  explicit  calculation  of  field-dependent 
interaction  integrals  with  straightforward  expressions  involving  the  residues  of  the  scat¬ 
tering  data,  the  method  provides  further  motivation  to  employ  inverse  scattering  methods 
in  the  design  of  optical  devices. 
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VII.  APPENDIX 

The  matrix  A(x)  : 

'0  1  0  0  0  O' 

0  0  f{rli)  a(c]  +  <%)  0  0 

0  0  0  0  /{■»)  n(cf  +  4) 

1  —x  e~mz  emz  e-*»r  e»&*  » 

0  ->  ie~n’  fa™ 

LO  0  fa-™  fa™  fa-™  faw  . 

where 

f{Vm)  =  (*?m  +  i&l)  (Vm  +  ifo)  (*?m  +  tfa)  m  =  1,2.  (93) 
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Table  1  Eigenvalues  and  corresponding  shape  factors  (eq.(90))  for  three  values  of  separation. 


Step  Index 

Truncated  Seek 2  ax 

High  -V  Profile 

2.26  x  109 

1.89  -7.56  x  106 

1.5  -  9.0  x  105 

Table  2  Representative  values  of  the  shape  factor  for  the  three  types  of  refractive  index  profile  considered  in 
this  paper.  The  step  index  profile  has  a  width  of  0.94  fim.  and  a  core  refractive  index  ni  =  3.1. 


List  of  Figures 

1.  Typical  multilayer  planar  waveguide  (two  layers  shown).  Guiding  regions  are  shown 

shaded. 

2.  Coupled  potentials  used  to  model  the  scattering  picture  of  weakly  coupled  planar 
waveguides. 

3.  Three  depressed-cladding  refractive  index  profiles  for  different  (01,02):  (0,0.001  a) 
(far  right),(0, 0.25a)  (dashed),  (0.499  a,  \/0.687  a)  (left). 

4.  Shape  factors  as  a  function  of  c\  for  C2  =  0.499  a  (top),  and  C2  =  0.25  a. 

5.  From  top  to  bottom:  the  area  under  the  potential,  the  shape  factor,  and  2  kp/i  as  a 
function  of  the  truncation  point. 
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Abstract.  Application  of  inverse  scattering  theory  for  designing  planar  optical  waveguides 
possessing  prescribed  propagation  constants  for  light  with  a  given  frequency  is  well  known. 
However,  waveguides  designed  using  such  a  method,  in  general,  will  not  be  able  to  transmit 
light  at  other  frequencies  with  the  same  propagation  constant.  In  order  to  overcome  this 
difficulty,  the  design  problem  for  te  modes  is  transformed  and  reformulated  to  an  equivalent 
inverse  problem  for  Schrodinger's  equation.  Then  using  inverse  scattering  theory,  the  potential 
as  a  function  of  a  modified  spatial  variable  is  recovered.  Next  the  important  problem  of 
finding  an  explicit  relation  between  the  actual  spatial  variable  and  the  modified  spatial 
variable  is  solved  and  a  systematic  procedure  is  developed  for  designing  waveguides  which 
have  the  same  propagation  constant  for  different  light  frequencies.  Existence  and  uniqueness 
questions  are  studied  and  some  model  calculations  are  presented. 


1.  Introduction 

Proper  values  of  propagation  constants  are  very  important  in  the  design  of  waveguides, 
since  they  govern  the  spatial  and  temporal  characteristics  of  the  signals  transmitted  in 
waveguides.  Systematic  proc  xiures  for  designing  waveguides  with  prescribed  propagation 
constants  appeal  to  the  existing  inverse  scattering  theories  [1-5],  which  were  originally 
developed  for  the  inverse  problems  in  quantum  mechanics.  In  standard  applications  of 
inverse  scattering  theories  for  designing  optical  waveguides  [6-9],  we  make  use  of  the  fact 
that  at  a  fixed  frequency  Maxwell’s  equations  governing  the  light  propagation  in  a 
waveguide  can  be  transformed  to  Schrodinger’s  equation  with  an  energy-independent 
potential.  In  this  equivalent  quantum  mechanical  inverse  problem,  the  bound  states 
energies  are  associated  with  the  prescribed  propagation  constants,  and  the  potential  is 
related  to  the  refractive  index  of  the  designed  waveguide. 

The  systematic  procedures  for  designing  waveguides  as  outlined  above  [6-9]  are 
applicable  as  long  as  we  are  interested  in  light  propagation  with  a  prescribed  propagation 
constant  at  a  single  frequency  through  the  designed  waveguide.  However,  such  a 
waveguide  in  general  will  not  have  the  designed  propagation  constants  for  light  with 
frequencies  other  than  the  specific  frequency  used  in  the  design  of  the  waveguide. 
Waveguides,  which  have  the  same  propagation  constants  for  different  light  frequencies, 
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have  important  applications  in  optics,  such  as  in  harmonic  generation,  wave  mixing, 
parametric  amplification,  and  multiplexing  [10-12]. 

The  need  for  developing  a  systematic  method  to  design  waveguides  having  the  same 
propagation  constant  for  different  light  frequencies  has  motivated  us,  in  this  preliminary 
study,  to  design  planar  optical  waveguides  that  have  the  same  propagation  constant  for 
te  modes  at  different  frequencies.  We  achieve  this  objective  by  showing  that  Maxwell’s 
equations  can  be  related  to  Schrodinger’s  equation  with  an  energy-dependent  potential 
and  that  the  requirement  for  the  waveguide  to  have  the  same  propagation  constant  for 
m  different  frequencies  is  shown  to  be  equivalent  to  the  corresponding  energy-dependent 
potential  supporting  m  bound  states  of  specified  values.  Having  reduced  the  problem  to 
an  inverse  Schrodinger  problem  for  energy-dependent  potentials,  we  then  subject  this 
Schrodinger  equation  to  a  transformation  [13-19]  which  reduces  the  inversion  to  a 
Schrodinger  inverse  problem  for  an  energy-independent  potential.  The  modified  inversion 
problem  is  then  solved  by  using  the  existing  Schrodinger  inverse  scattering  methods  in 
one  dimension  [1-5].  However,  since  the  equivalent  inverse  Schrodinger  problem  is 
formulated  with  respect  to  a  modified  spatial  variable  and  not  the  actual  spatial  variable, 
the  energy  independent  potential  found  will  not  be  of  any  use  unless  the  connection 
between  the  actual  spatial  variable  and  the  modified  spatial  variable  is  established.  We 
study  this  important  problem  in  detail  and  find  an  explicit  relation  between  the  actual 
and  the  modified  variable,  which  then  enables  us  to  make  use  of  the  energy  independent 
potential  and  develop  a  systematic  and  practical  procedure  to  design  waveguides  which 
have  the  same  propagation  constant  for  different  light  frequencies. 

In  section  2  we  review  the  problem  of  electromagnetic  wave  propagation  in  a  planar 
waveguide.  Section  3  deals  with  transforming  Maxwell’s  equations  to  Schrodinger’s 
equation  with  an  energy-independent  potential  and  developing  a  systematic  method  to 
design  a  waveguide  which  has  the  same  propagation  constant  for  different  frequencies. 
In  section  4  examples  of  practical  interest  in  waveguide  design  are  presented.  We  find 
that  the  proposed  method  enables  us  to  design  waveguides  which  have  the  same 
propagation  constant  for  any  finite  number  of  different  light  frequencies.  The  procedure 
leads  to  solutions  which  depend  on  infinitely  many  arbitrary  parameters.  Of  course,  this 
non-uniqueness  can  be  further  manipulated,  enabling  the  designed  waveguide  to  have 
other  desirable  properties. 


2.  Statement  of  the  problem 

Propagation  of  electromagnetic  waves  in  a  planar  optical  waveguide,  with  refractive 
index  varying  continuously  only  in  one  direction  say  x,  is  analyzed  by  assuming  that  the 
electric,  E,  and  magnetic,  //,  fields  have  the  following  forms  [6]: 

E,(x,y,z,t)  =  Em(xyd*"~f*  (2.1) 

(2.2) 

where  x,y,  and  i  are  the  cartesian  coordinates  with  z  along  the  axis  of  the  waveguide, 
t  is  time,  a  represents  the  components  of  a  vector  in  the  x,y,  or  z  directions,  c o  is  the  light 
frequency,  and  /?  is  the  propagation  constant.  Subsitution  of  (2. 1 )  and  (2.2)  in  Maxwell’s 
equations  lead  to  the  following  equation  [6], 

\l/(x)  +  [nJ(x,  k0)kl  -  P^ix)  =  0  (2.3) 
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for  the  te  modes.  In  (2.3)  the  positive  function  n(x,k0)  is  the  refractive  index,  k0  =  co/c , 
with  c  being  the  speed  of  light  in  vacuum,  and  if/(x)  is  the  field  function  associated  with 
the  electromagnetic  fields  under  consideration.  The  field  function  tp(x)  is  to  decay  fast 
enough,  as  |x|  increases,  so  that  the  field  function  is  associated  with  finite  energy  which 
is  mostly  confined  to  the  inside  of  the  waveguide. 

It  is  well  known  that  the  differential  equatio  (2.3)  with  the  above  condition  can  have 
solution  only  for  certain  values  of  /?,  which  are  called  the  eigenvalues  of  the  differential 
equation,  and  for  the  problem  at  hand  correspond  to  different  possible  propagation 
constants  of  the  waveguide.  Therefore  in  the  waveguide  design  problem  at  a  fixed 
frequency,  one  only  need  to  find  n(x,  k0)  which  is  associated  with  the  desired  propagation 
constants  P  for  the  specified  frequency.  A  standard  procedure  to  solve  this  design 
problem  is  to  appeal  to  the  theory  of  inverse  scattering  which  was  first  developed  in 
quantum  mechanics,  where  one  has  to  find  the  potential  from  the  spectral  data  of  the 
associated  Schrodinger  equation  [1-5].  In  order  to  be  able  to  make  use  of  the  well 
developed  methods  of  inverse  scattering  theory  in  quantum  mechanics  [1-5],  one 
transforms  (2.3)  into  a  Schrodinger  differential  equation  form 

^3  ip(x)  +  [k2  -  k2  V(x,  k„  )]4>(x)  =  0  (2.4) 

where 

k2  =  ni(k0)kl  -  P2  (2.5) 

V(x,  k0)  =  ni (kQ)  -  n2(x,  k0)  (2.6) 

with  nx(k0)  being  the  refraction  index  for  |x|  -♦  oo. 

Having  transformed  the  Helmholtz  equation  (2.3)  into  a  Schrodinger  equation  (2.4), 
one  notes  that  the  design  problem  of  optical  waveguide,  that  is  finding  the  index  of 
refraction  n(x,  k0)  which  gives  us  the  desired  propagation  constant  P  for  the  specified  k0, 
is  reduced  to  an  inverse  scattering  problem  in  quantum  mechanics,  where  the  potential 
k\  V{x,k0)  is  to  be  deduced  from  the  information  on  the  bound  states  and  the  reflection 
coefficients.  In  this  quantum  mechanical  formulation  of  the  problem,  one  refers  to  k2  as 
the  energy  of  the  system  and  the  eigenvalues  as  the  bound  state  energies,  which  we  will 
denote  by  —  y2 .  Of  course  as  can  be  seen  from  (2.5)  these  binding  energies,  y2 ,  are  related 
to  the  desired  propagation  constants  through  the  relation 

r  =  p2  -  nl(k0)k2.  (2.7) 

From  (2.7)  it  follows  that  specification  of  the  propagation  constant  p  and  frequency 
tu,  will  give  us  the  needed  bound  state  energy  information  for  the  analogue  quantum 
mechanical  problem.  Having  established  the  connection  between  the  optical  waveguide 
and  the  inverse  quantum  mechanical  problem,  it  is  then  straightforward  to  use  existing 
methods  [5-7]  to  find  the  desired  potential  kl  V(x,  k0)  associated  with  the  bound  states 
and  the  reflection  coefficients  and  then  find  the  required  refractive  index  n(x,k0)  from 
F(x,  fc0)  using  (2.6).  However,  this  standard  approach  is  useful  for  designing  waveguides 
associated  with  only  one  fixed  frequency.  That  is,  since  the  inversion  potential  k%  V(x ,  k0) 
depends  on  the  frequency  w,  if  we  change  c o,  the  potential  will  change  resulting  in  change 
of  the  binding  energy  y2 .  In  other  words  the  waveguide  designed  will  not  have  the  desired 
propagation  constant  P  at  other  frequencies.  Therefore  if  we  are  interested  in  designing 
waveguides  which  have  the  same  propagation  constant  for  different  frequencies,  the 
method  as  stated  above  is  not  able  to  provide  us  with  the  desired  profile.  We  will  show 
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in  the  next  section  that  it  is  still  possible  to  use  the  results  of  inverse  scattering  theory 
[1-5]  to  design  waveguides  which  can  have  the  same  propagation  constant  for  different 
light  frequencies. 

Before  concluding  this  section,  let  use  note  that  the  refractive  index  n(x,l fco)  in  general 
is  a  function  of  both  the  spatial  variable  x  and  also  the  wavenumber  k0.  The  dependence 
of  the  refractive  index  on  k0  or  the  frequency  of  the  light  propagating  through  the 
waveguide  is  a  very  interesting  and  important  topic.  However,  in  this  preliminary  study, 
for  the  sake  of  simplicity  in  presentation,  we  restrict  the  study  to  refractive  indexes  which 
are  twice  differentiable  with  respect  to  x  and  have  the  following  type  of  dependence  on 
k0  and  x: 

n(x,k0)  =  na,(k<))ti(x)  (2.8) 

where  tj(x)  is  a  function  of  x  only  and,  which  tends  to  1  as  |x|  tends  to  infinity.  Since 
noo(ko)  *s  associated  with  the  refractive  index  of  the  cladding,  it  will  be  assumed  that 
nK(k0)  is  an  arbitrary  but  known  function  of  k^.  In  other  words,  in  this  study  we  make 
the  assumption  that  the  refractive  index  is  made  up  of  a  known  positive  function  nT  (fc0), 
multiplied  by  a  positive  function  rj(x)  which  is  a  function  of  x  only.  In  this  study  we  also 
need  to  restrict  rj(x)  to  class  of  function  which  satisfy  the  following  inequality: 

f  toM-  l|d*<  oo.  (2.9) 

The  design  problem  to  be  presented  in  section  3  is  to  develop  a  systematic  procedure  for 
finding  the  frequency-independent  function  ri(x)  corresponding  to  a  refractive  index 
which  will  allow  different  light  frequencies  to  propagate  through  the  waveguide  with  the 
same  propagation  constant  p. 


3.  The  inversion  procedure 

As  was  shown  in  the  previous  section,  our  design  problem  is  to  develop  a  systematic 
procedure  to  design  waveguides  which  have  the  same  propagation  constant  p  for  all 
different  light  frequencies  of  interest.  In  order  to  be  able  to  restate  this  design  specification 
in  the  equivalent  inverse  quantum  mechanical  problem  in  a  more  transparent  fashion, 
let  us  rewrite  the  differential  equation  (2.4)  in  the  following  manner: 

—  *(x)  +  [k2  -  k1  K,(x)  -  V2{x)W(x)  =  0  (3.1) 

where 

Vi (*)  =  k0)/n2M  =  1  -  n2W  (3.2) 

V2(x)  =  p2Y,(x).  (3.3) 

It  should  be  remembered  that  throughout  the  discussion,  the  propagation  constant 
P  is  fixed  but  the  frequency  cu  can  take  different  values,  (o,  with  i  =  l ,2 . . .  ,m.  Since  (3.1) 
is  the  same  as  (2.3),  the  eigenvalues  of  (3.1)  will  be  the  same  as  those  of  (2.3)  and  will 
be  related  to  frequency  co  and  propagation  constant  P  through  (2.5).  However,  the 
advantage  of  writing  (2.3)  in  the  form  of  (3.1)  is  the  fact  that  (3.1)  clearly  shows  that  our 
design  problem  is  equivalent  to  an  energy-dependent  Schrodinger  inverse  problem  and 
we  are  interested  in  finding  V;(x)  and  V2(x)  when  binding  energies  of  (3.1)  are  specified 
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according  to  following  equation: 

y-=P2-  nl (*,*)*£  /  =  1 ,2, . . . ,  m  (3.4) 

where  k0,  =  a),/c.  Equation  (3.4)  is  nothing  but  (2.7),  and  it  is  only  written  to  emphasize 
the  fact  that  in  the  problem  of  interest  we  are  not  given  a  single  bound  state  energy,  as 
may  be  inadvertently  deduced  from  the  fact  that  we  only  have  one  propagation  constant, 
but  in  fact  we  are  given  m  bound  states  energies.  These  binding  energies  can  be  easily 
computed  from  (3.4)  by  substituting  the  desired  different  values  of  frequency  to, 
which  we  would  like  to  propagate  through  the  waveguide  with  the  same  propagation 
constant  /?. 

Equation  (3.1)  as  it  stands  is  not  in  the  standard  Schrodinger  equation  form  and 
therefore  existing  inversion  methods  for  energy-independent  potentials  cannot  be  directly 
applied.  However,  similar  equations  have  been  dealt  with  when  one  tries  to  solve  inverse 
problems  for  angular-momentum-dependent  potentials  [13-14]  and  wave  equations  in 
one  dimension  [15-19].  Motivated  by  these  results,  let  us  then  transform  our  independent 
variable  x  to  p,  through  the  following  relation: 

P(x)  =  J '  dr  Vl  -  V,(t)  =  £  n(t) At.  (3.5) 

In  view  of  the  fact  that  rj(x)  is  a  positive  function,  the  above-defined  mapping  is 
one-to-one  and  the  inverse  mapping  exists.  This  allows  us  to  write  the  quantities  of 
interest  as  either  functions  of  x  or  as  functions  of  p,  depending  on  which  representation 
is  more  suitable  for  solving  the  inversion  problem.  With  this  observation  in  mind  let  us 
define  a  modified  field  function  through  the  relation 

<f>(x)  =  Ip(x)l<x(x)  (3.6) 

with 

=  yjlhix)-  (3.7) 

Changing  variable  in  (3.1)  from  x  to  p  and  making  use  of  (3.6)  one  can  rewrite  (3.1) 
in  the  following  form: 

^2  ^(P)  +  1*2-  W(P))$(P)  =  0  (3.8) 


where 


W(p)  =  [^s  ^)][2^(p)]-'  -  i/(p)J  [2 i}(p)]-2  -  P7[  1  -  1  InipH  (3.9) 

In  (3.9)  $(p)  and  rj(p)  are  the  field  function  <p(x)  and  the  refractive  index  tj(x),  written 
as  functions  of  p,  respectively. 

The  advantage  of  working  with  (3.8)  instead  of  (3.1)  is  clear.  Equation  (3.8)  is  the 
Schrodinger  equation  for  the  energy-independent  potentials,  whose  inverse  problem  is 
well  studied.  Furthermore,  let  us  note  that  if  $  is  an  eigenfunction  of  (3.8),  then  the 
associated  field  function  ip  is  also  an  eigenfunction  of  (3. 1 ).  In  other  words  the  eigen¬ 
values  of  the  two  equations  are  the  same  and  therefore  the  design  problem  is  reduced  to 
finding  W(p)  with  bound  state  energies  specified  by  (3.4).  This  is  the  classical  inverse 
quantum  mechanical  problem  and  the  solution  to  it  is  well  known  [1-6]-  The  only  point 
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that  we  need  to  emphasize  is  that  the  solution  is  not  unique  and  even  if  we  specify  not 
only  the  bound  state  energies  but  also  the  reflection  coefficients,  still  the  inversion  result 
we  depend  on  m  arbitrary  parameters,  which  in  the  design  problem  could  be  used  to  our 
advantage.  However,  for  the  moment  let  us  assume  that  a  potential  W(p)  associated  with 
the  given  bound  states  has  been  obtained  and  the  refractive  index  rj(p)  solution  of  the 
nonlinear  differential  equation  (3.9)  with  the  boundary  conditions 

lim  fj(p)  =  1  (3.10) 

|/H  -*oo 

has  been  found  as  a  function  of  the  intermediate  variable  p.  Then  the  only  remaining 
problem  is  to  find  the  refractive  index  as  a  function  of  the  spatial  variable  x,  when  the 
refractive  index  as  a  function  of  p  is  known.  In  order  to  achieve  this  objective  we  make 
use  of  (3.5)  to  deduce  the  following  relation: 


Now  since  fj(p)  is  a  known  function,  equation  (3.1 1)  can  be  used  to  find  x  as  function 
of  p.  In  other  words  the  one-to-one  function  F{p)  can  be  computed  and  its  inverse  F' 1  (x) 
can  also  be  found.  Noting  that  p  =  F~'(x),  we  are  then  in  a  position  to  find  tj(x)  — 
ij(F~'(x)).  By  construction  the  so-designed  waveguide  will  have  the  same  propagation 
constant  P  for  all  light  frequencies  to,,  with  i  =  1,2, ...  ,m. 

In  principle  the  above  procedure  enables  us  to  design  waveguides  with  the  same 
propagation  constants  P  for  different  frequencies  provided  that  we  can  find  ij(p).  In  order 
to  show  the  existence  of  the  solution  to  (3.9)  and  develop  a  practical  method  for  finding 
the  solutions,  we  note  that  Berryman  and  Greene  [18],  in  dealing  with  inverse  problems 
for  elastic  waves,  have  shown  that  the  impedance  can  be  either  recovered  directly  by 
solving  a  linear  second-order  differential  equation,  which  can  be  regarded  as  the 
analogue  of  (3.9),  or  indirectly  by  working  with  the  wavefunction  associated  with  zero 
frequency.  Motivated  by  this  result  [18],  let  us  study  the  wavefunctions,  solutions  to  (3.1) 
and  (3.8),  at  zero  frequency  which  corresponds  to  k?  —  —  ft2.  We  note  that  when  tw  =  0 
equation  (3.1)  simplifies  to  the  following  equation: 

^(*)-/^(*)  =  0  (3.12) 

with  e±fx  being  its  two  fundamental  solutions.  Let  us  also  denote  <j>±  (p)  as  the  solutions 
to  (3.8)  for  k2  =  —  p2  with  asymptotic  behaviours  of  the  form  cifp  for  p  tending  to  +  oo, 
respectively.  We  should  note  that  $  t  (p)  are  linearly  independent.  Otherwise,  we  are 
forced  to  accept  that  —  P2  is  an  eigenvalue  of  (3.8).  However,  this  is  not  the  case  since 
we  are  assuming  that  the  potential  W(p)  is  chosen  in  such  a  way  that  (3.8)  has  eigenvalues 
—  y2  as  given  by  (3.9).  Also  we  assume  W(p)  is  such  that  the  solution  $  ±  (p)  to  (3.8)  exist 
for  all  real  values  of  p,  and  the  associated  function  q(x)  satisfies  (2.9).  Having  defined  the 
desired  solutions  to  (3.1)  and  (3.8)  for  k2  =  -  p2 ,  we  then  make  use  of  the  relation  (3.6) 
to  find 

i±<P)-y/W)Atc*''  (3.13) 

where  At  =  exp(±/Jjo  00  M-x)  -  l]dx).  It  should  be  noticed  that  (3.13)  is  the  analogue  of 
equation  (40)  of  Berryman  and  Greene  { 1 8],  however,  in  order  for  (3. 1 3)  to  be  of  practical 
use  we  need  to  eliminate  its  dependence  on  x  by  taking  the  derivative  of  (3.13)  with 
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respect  to  p.  Making  use  of  the  existing  relation  between  p  and  x  a»  given  by  (3.1 1),  we 
obtain 


^h(p)  =  [^(P)  ±  2py±(p)[2i(p)]-'  (3.14) 

which  can  be  viewed  as  two  linear  first-order  differential  equations  satisfied  by  rjipi  The 
solution  to  (3.14)  can  be  written  as 

n(p)  =  2p$l(p)  r«?;J(r)d/  (3.15) 

Jp 

or 

n(p)  =  2p$l(p)  J"^:2(/)dr.  (3.16) 

It  should  be  noted  that  if  ij(p)  defined  by  (3. 1 5)  becomes  different  from  fj(p)  as  defined 
by  (3.16),  then  from  equations  (3.13)  and  (3.14)  one  is  forced  to  accept  that  the  linear 
second-order  differential  equation  (3.8)  has  more  than  two  linearly  independent 
solutions.  Of  course  this  not  being  the  case,  proves  that  the  refractive  index  tj(p)  as 
defined  by  (3. 1 5)  or  (3.16)  are  identical  and  one  can  use  either  representation  to  compute 
the  refractive  index.  Also,  since  we  are  finding  the  refractive  index  in  such  a  roundabout 
way,  one  is  justified  in  asking  whether  this  refractive  index  actually  satisfies  (3.9). 
Performing  the  necessary  operations,  it  is  very  easy  to  verify  that  indeed  rj(p)  as  defined 
by  (3.15)  or  (3.16)  satisfies  (3.9).  Furthermore,  using  a  similar  proof  to  that  developed 
to  show  (3.15)  and  (3.16)  lead  to  the  same  refractive  index,  we  can  conclude  that  the 
solution  to  (3.9),  satisfying  the  boundary  conditions  specified  by  (3.10),  is  unique  and  is 
given  by  (3.15)  or  (316). 

Positivity  of  fj(p)  can  be  deducted  from  (3.13),  (3.15)  or  (3.16).  Equations  (3.15)  and 
(3.16)  show  that  tj(p)  is  non-negative.  To  show  positivity  of  let  us  assume  that  there 
exist  a  point  p0  such  that  fj(Po)  =  0.  Then  (3.13)  will  imply  that  <f>±(p0)  =  0.  From  this 
information  we  can  deduce  that  the  Wronskian  of  $+ (p)  and  $  _  (p)  is  equal  to  zero.  In 
other  words  the  solutions  $±  (p)  of  (3.8)  are  linearly  dependent.  Then  we  appeal  to  the 
fact  that  $±(p)  are  linearly  independent  and  conclude  that  the  point  p0  such  that 
i)(p0)  =  0  does  not  exist  and  i\(p)  is  a  positive  function. 

The  proposed  nr:thod  for  finding  the  refractive  index  rj(p),  seems  to  have  replaced  the 
need  for  finding  the  solution  to  the  nonlinear  boundary  value  problem  as  given  by  (3.9) 
by  the  need  to  find  the  solution  to  the  linear  equation  (3.8)  for  a  special  value  of  the 
energy  k 2 .  Although  this  by  itself  is  a  great  simplification,  it  should  be  noted  that  the  gain 
is  even  greater  when  we  remember  that  any  standard  inversion  procedure  which  we  use 
to  find  W(p)  will  also  be  able  to  give  as  the  wavefunction  for  different  k2  values  without 
having  to  directly  solve  the  associated  Schrodinger  equation.  In  other  words  the  functions 
$±  (p)  needed  for  calculation  of  rj(p)  can  be  easily  found  and  we  will  not  need  to  appeal 
to  numerical  methods  to  solve  the  linear  differential  equation  (3.8)  for  k2  ~  —  fi2 . 
Furthermore,  the  solution  given  in  the  form  of  (3.15)  and  (3.16)  enable  us  to  easily 
integrate  (3.11)  and  find  the  dependence  of  the  spatial  variable  x  on  the  intermediate 
variable  p: 
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or 


<3is) 

Again  the  two  representations  of  x  as  given  by  (3.17)  and  (3. 18)  are  equivalent,  and  can 
be  used  to  relate  the  spatial  variable  x  to  the  intermediate  variable  p. 

We  conclude  this  section  by  noting  that  the  proposed  procedure  is  able  to  give  us  a 
practical  method  for  designing  waveguides  which  have  the  same  propagation  constant 
for  different  light  frequencies.  The  procedure  does  not  lead  to  a  unique  solution  and  this 
of  course  is  of  practical  importance  since  it  enables  us  to  design  waveguides  having 
further  desirable  properties.  The  sources  of  non-uniqueness  are  due  to  the  fact  that  in 
this  design  problem  the  reflection  coefficients  are  not  specified  and  can  be  chosen 
arbitrarily  and,  furthermore,  for  each  required  bound  state  we  also  have  a  normalization 
parameter  which  is  arbitrary.  In  order  to  illustrate  the  procedure  in  more  detail  and  see 
some  of  the  effects  of  the  existing  arbitrariness  in  the  procedure,  in  section  4  we  present 
examples  which  are  also  of  practical  interest  in  waveguide  design. 


4.  Examples 


In  this  section  we  study  the  design  of  a  waveguide  which  allows  two  frequencies  cw,  and 

(1)2  to  propagate  with  the  same  propagation  constant  fi.  Following  the  procedure 

developed  in  section  3  we  first  use  (3.4)  to  define  the  bound  state  energies  associated  with 
this  problem. 

yf  =  -  "l>(fcoi)*2i  (4-i) 

y]  =  F-  nl(k01)k2O2  (4.2) 


where  n^(km)  and  n^fa 2)  are  the  refractive  index  of  the  cladding  at  the  frequencies  co, 
and  o)2,  and  k0l  —  /c  and  k^  =  o)2/c.  Having  defined  the  desired  bound  state  energies, 
we  are  now  ready  to  appeal  to  the  well  known  results  of  "ay  and  Moses  [20]  to  find  the 
bound  state  wavefunctions  $,(p)  and  $2(p)  and  the  associated  potential  W(p): 


and 


$i(p) 

$2(p) 


Axt"'  AyA2{yy-y2^^ 
A(p)  2y2(y,  -1-  y2)A(p) 

A2cw  A,A2(  y,-y2)cl*'  +  ^ 

A(p)  2y,(y,  +  y2)A(/>) 


W(p)  =  2 


d_ 

dp 


+  $2(p)<:TlP\ 


(4.3) 

(4.4) 

(4.5) 


where 


i4,e 


A2e**  .  AyA2(yy  -  y2)Je 


,2  2(7 1  +  7;V> 


A(p)  =  I  +  - - 1-  - -  h  - - - - - - - Tf 

2y,  2y2  4  yty2(y,  +  y2r 

and  A,  and  A2  are  arbitrary  positive  constants. 


(4.6) 
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fco |  =  2.0270 //m'1  and  ^  =  4.0537 /an~‘. 


Having  found  the  potential  and  the  bound  state  wavefunctions  we  are  then  in  a 
position  to  use  equation  (2.1)  of  Kay  and  Moses  [20]  to  find  the  needed  wavefunctions 
<f>±(p)  for  k?  =  —p1  without  having  to  solve  the  differential  equation  (3.8)  directly: 


and 


$  +  <J>)  = 

#-(p)  = 


3 1  W 

V.  +  0 


&(p)eT>1’  la. 

7:  +  /*  J 


fa  ~  P  J 


(4.7) 

(4-8) 


It’s  easy  to  verify  that  the  so-defined  <j>t(p)  have  the  desired  asymptotic  behaviours  and 
are  solutions  to  (3.8).  Therefore  they  can  be  used  in  (3.15)  or  (3.16)  and  (3.17)  or  (3.18) 
to  find  the  refractive  index  and  the  spatial  variable  x  as  a  function  of  the  intermediate 
varaible  p  numerically.  Having  found  and  x  =  F\p),  the  procedure  is  then  complete 
and  the  refractive  index  ij(jc)  can  be  numerically  obtained.  The  result  of  the  numerical 
computations  are  presented  for  different  values  of  A,,  A2  and  different  light  frequencies 
in  figures  1-4.  In  these  examples  we  have  assumed  that  the  refractive  index  n(x,  k^)  as 
defined  by  (2.8)  is  independent  of  wavenumber  k0  and  can  be  written  as  t](x)nx. 

It  should  be  emphasized  that  the  examples  presented  here  are  only  for  the  sake  of 
demonstrating  the  proposed  method.  Practical  implementation  of  this  technique  and 
actual  fabrication  of  such  waveguides  need  further  study.  Also,  for  the  sake  of  simplicity 
in  presentation,  we  have  only  used  reflectionless  potentials  in  these  examples.  However, 


«(^m] 


Figure  2.  Graph  of  the  refractive  index 
rt(x,k0 )  =  (f(.xVx  as  defined  by  (2.8)  and 
(3.15).  Symbols  and  the  constants  are  the 
same  as  m  figure  I . 
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PlMm] 


Figure  3.  Plot  of  fl  =  x  -  p  = 
F[p)  -  p.  The  solid  line  corresponds 
to  A,  =  Aj  =  2  and  the  dotted  line 
is  associated  with  ^,=4  and 
A1  =  8.  In  all  cases  p  =  6.061 1  pm~\ 
«.(*<>)  =  «.  =  1  4850,  *,,= 

’i&lTQpm''  and  =  4.0537 /on-1. 


in  actual  applications  of  the  method,  we  should  remember  that  any  potential  which  has 
the  proper  bound  state  energies,  including  those  which  are  not  reflectionless,  can  be  used. 
Such  potentials  can  be  found  by  appealing  to  the  Faddeev-Marchenko  method  [5,21]. 

W{p)  =  2±-K{p,p)  (4.9) 


where 


K(p  +  Q  +  M(p,Q+  fP  K(p,Om  +  0^  =  0  C  <P  (4.10) 

J  -  00 

and 

M(p)  =  (l/27t)  f  "  d k  R(k)e-'k‘’  +  £  At  ew  (4.11) 

with  R(k)  being  the  reflection  coefficient.  Having  found  W(p),  one  can  apply  the 
proposed  method  to  find  the  refractive  index  associated  with  potentials  which  are  not 
reflectionless.  Use  of  potentials  with  R(k)  #  0  may  be  preferable,  if  such  potentials  lead 
to  waveguides  with  refractive  indexes  which  are  easier  to  fabricate. 

The  above  equations  show  that  in  order  to  find  W(p)  uniquely,  we  not  only  need  the 
bound  state  information  and  the  normalization  constants  A„  but  also  we  need  to  know 
the  reflection  coefficients  R(k)  for  all  real  values  of  k.  In  view  of  the  fact  that  in  fibre 
optics  design  usually  only  the  value  of  propagation  constants  are  specified  and  R(k)  is 


Figure  4.  Graph  of  the  refractive 
index  n(x,k^)  *  ti(x)nm.  Symbols 
and  the  constants  are  the  same  as  in 
figure  3. 
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not  given,  it  then  follows  that  such  a  design  problem  involves  a  great  degree  of  non¬ 
uniqueness.  A  very  interesting  question  that  was  raised  by  the  referees  is  with  regard  to 
the  nature  of  this  non-uniqueness  and  its  physical  implications.  Let  us  study  this  question 
by  noting  that  in  order  to  remove  this  non-uniqueness,  we  need  to  know  R(k)  for  all  real 
values  k.  In  view  of  the  definition  of  k  as  given  by  (2.5),  we  note  that  the  wavevectors 
of  interest,  k  —  (kx,  k,),  will  fall  into  one  of  the  following  two  categories.  Case  (a)  is  when 
/?  is  real  and  n1 2 3 4 5 6 7 8 9 10x(k0)kl  >  /?2,  resulting  in  both  components  of  At  to  be  real.  That  is 
At,  =  k  =  +  [«i(At0)Ato  —  P1]'11  and  Ac.  =  /?.  Case  (b)  corresponds  to  the  evanescent  waves 
[22],  where  /?  =  Ac.  =  —  i/5  is  purely  imaginary  but  kx  =  k=  ±  [r4(Ac0)fcj;  +  j32],/2  is  still 
real.  From  the  above  analysis,  it  also  follows  that  in  case  (a)  for  large  values  of  |x|  the 
wave  will  behave  like  a  free  wave  and  therefore  from  the  point  of  view  of  geometric  optics 
it  would  correspond  to  refracted  rays  {22],  This  analysis  shows  that  data  on  R(k)  are 
associated  with  waves  which  are  significant  only  in  the  spatial  transient  region  and  their 
powers  are  significantly  diminished  in  the  spatial  steady-state  region  [22).  The  only  waves 
that  will  have  significant  power  for  large  values  of  z,  that  is  in  the  spatial  steady-state 
region,  are  the  bound  waves  (22).  Of  course  propagation  constants  of  such  waves,  /?,, 
have  played  a  very  important  role  in  our  design  procedure.  This  observation  clarifies  the 
nature  of  the  existing  non-uniqueness  in  our  design  problem.  It  shows  that  the  main 
difference  between  the  different  waveguides  which  can  be  deduced  from  the  proposed 
method  is  in  their  radiation  properties  in  the  spatial  transient  region,  which  is  usually  a 
short  distance  from  the  source.  However,  for  most  of  the  length  of  the  proposed 
waveguides,  that  is  in  the  spatial  steady-state  region,  waves  associated  with  data  R(k)  will 
not  be  significant  and  only  the  bound  waves  will  be  present.  In  other  words,  in  the  spatial 
steady-state  region,  all  of  the  proposed  waveguides  will  perform  similarly  as  far  as  the 
bound  waves  are  concerned.  It  should  again  be  emphasized  that  although  we  are  mainly 
interested  in  the  propagation  of  the  bound  waves,  the  existing  non-uniqueness  can  play 
an  important  role;  such  as  the  ease  of  fabrication  of  the  waveguide  or  coupling  of  energy 
from  the  source  to  the  waveguide.  Of  course  such  a  study  is  beyond  the  scope  of  this 
paper  but  it  deserves  further  consideration  both  for  planar  waveguides  and  circular 
waveguides  [6],  where  the  same  type  of  non-uniqueness  also  exists. 
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We  have  formulated  a  matrix  eigenvalue  problem  for  cylindrical  optical  libers  from  a  set  of  finite  difference  equations.  Numer¬ 
ical  solution  of  this  problem  yields  the  propagation  constants  for  propagating  modes.  The  method  can  be  used  for  arbitrary  index 
profiles,  does  not  require  the  explicit  evaluation  of  Bessel  or  modified  Bessel  functions,  and  does  not  use  iterative  methods  to 
search  for  the  propagation  constants  as  was  the  case  in  earlier  proposed  methods  using  finite  differences.  The  method  is  accurate, 
fast,  and  simple.  We  have  established  the  convergence  and  stability  of  this  method,  and  explored  the  effects  of  finite  cladding 
width  on  the  dispersion  characteristics. 


1.  Introduction 

Wave  propagation  in  optical  fibers  has  been  ana¬ 
lyzed  using  various  methods.  We  will  be  using  a  fi¬ 
nite  difference  method.  Other  methods  proposed  to 
find  the  propagation  constants  of  guided  modes  in 
optical  fibers  with  arbitrary  refractive  index  profiles 
include  the  WKBJ  method,  variational  method, 
power  series  expansion  method,  staircase  approxi¬ 
mation  method,  and  finite  element  method. 

The  WKBJ  method  (1,2]  is  a  geometrical  optics 
approximation  that  works  whenever  the  refractive 
index  of  fiber  varies  only  slightly  over  distances  of 
the  order  of  the  optical  wavelength  and  are  appli¬ 
cable  only  to  thick  fibers  in  which  many  modes  can 
propagate.  For  those  fibers  in  which  only  a  few  modes 
propagate,  the  error  of  the  WKBJ  method  increases 
intolerably  and  this  method  is  not  applicable  to 
modes  near  cutoff.  Besides,  the  effect  of  an  index 
valley  at  the  core-cladding  boundary,  which  plays  an 
important  role  in  reducing  multimode  dispersion, 
cannot  be  treated  by  the  ordinary  WKBJ  method.  In 


the  variational  method  the  scalar  wave  equation  is 
converted  into  a  variational  problem  subject  to  the 
given  boundary  conditions.  The  variational  problem 
is  solved  either  by  using  the  Rayleigh-Ritz  method 
(3)  or  perturbation  method  [4],  In  the  Rayleigh- 
Ritz  method  the  eigenfunction  is  expressed  in  terms 
of  a  set  of  orthogonal  functions  and  the  variational 
function  is  minimized.  The  disadvantage  is  that  we 
need  to  assume  a  trial  function  ( 5],  In  the  pertur¬ 
bation  method  of  analysis,  the  computation  of  the 
propagation  characteristics  for  an  arbitrary  profile  is 
done  by  correcting  the  solution  for  a  uniform  core 
fiber  considering  the  difference  in  the  profile  as  the 
perturbation  term. 

The  power  series  expansion  method  [6]  consists 
of  expressing  the  refractive  index  for  the  field  term 
by  term.  This  method  is  useful  only  for  cases  in  which 
the  refractive  index  profile  can  be  expressed  by  a  rel¬ 
atively  simple  power  series.  In  some  cases  the  series 
do  not  converge  and  this  method  is  not  applicable 
[7].  In  the  staircase  approximation  method  [8,9] 
the  refractive  index  is  approximated  by  an  appro- 
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priate  staircase  function.  The  wave  equation  is  solved 
in  each  stratified  layer  and  the  solutions  are  then 
connected  at  the  cylindrical  boundaries  between  these 
layers  to  obtain  the  proper  solution  representing  the 
propagation  characteristics.  The  number  of  layers 
should  be  infinite  in  order  that  the  refractive  index 
profile  approaches  that  of  actual  fiber  profiles.  Thus, 
the  results  of  the  propagation  constant  will  differ  from 
the  actual  values  when  a  finite  number  of  layers  is 
used.  A  large  number  of  layers  requires  considerable 
computer  time  and  hence  in  this  method  the  accu¬ 
racy  and  computer  time  are  traded  off. 

The  fiber  problem  has  been  analyzed  by  Okamato 
and  Okoshi  [10]  using  a  finite  element  method  for¬ 
mulated  in  the  axial  fields.  The  problem  with  this 
method  is  that  it  suffers  from  spurious  modes  when 
the  finite  elements  are  not  carefully  chosen  [11]. 
Lenahan  [12]  has  formulated  a  matrix  eigenvalue 
problem  from  a  finite  element  analysis  using  the 
Galerkin  weighted  residual  method.  To  achieve 
computational  efficiency,  a  piecewise  linear  approx¬ 
imation  to  the  solution  function  must  be  used. 

In  this  paper  we  present  an  efficient  finite  differ¬ 
ence  method  to  find  the  propagation  constants  of  op¬ 
tical  fibers  with  arbitrary  refractive  index  profiles. 
The  method  does  not  involve  a  search  procedure  to 
find  the  propagation  constants,  nor  does  it  require 
explicitly  evaluating  Bessel  and  modified  Bessel 
functions,  as  was  the  case  in  the  earlier  works  on  fi¬ 
nite  difference  analysis  of  optical  fibers  [13,14].  We 
construct  a  matrix  equation  from  a  set  of  simulta¬ 
neous  finite  difference  equations  governing  the 
propagation  in  an  optical  fiber  and  solve  for  the  ei¬ 
genvalues  to  obtain  the  propagation  constants.  In 
sect.  2  we  give  the  mathematical  formulation  of  the 
discretized  differential  equation  at  various  grid  points 
in  the  radial  direction  and  the  construction  of  a  ma¬ 
trix  equation  incorporating  the  boundary  conditions 
at  the  core-cladding  interface  and  the  jacket  Ex¬ 
tending  our  method,  which  is  formulated  for  a-in- 
dex  fibers,  to  arbitrary  refractive  index  profiles  is 
covered  in  sect.  3.  In  sect.  4  we  discuss  the  numerical 
evaluation  of  propagation  constants  and  present  re¬ 
sults.  This  includes  a  discussion  of  the  convergence 
and  stability  of  the  method  along  with  the  effect  of 
the  number  of  grid  points  on  the  computation,  and 
the  effects  of  finite  cladding  width  on  dispersion 

394 


characteristics.  Our  conclusions  are  given  in  sect.  S. 


2.  Mathematical  formulation 

The  optical  fibers  considered  here  are  inhomoge¬ 
neous  dielectric  cylinders  of  radius  a  called  the  "core” 
surrounded  by  a  homogeneous  refractive  index  me¬ 
dium  called  the  “cladding”.  The  cladding,  in  turn,  is 
encased  in  a  highly  lossy  material  called  the  jacket. 
A  representative  fiber  cross-section  is  shown  in  fig.  1. 

The  refractive  index  profile  of  the  fiber,  called  an 
a-index  profile,  is  given  by 

n2(r)=n\[\-2pA(r/a)a\  ,  forO <r<a, 

=  n,[l-2i4],  forr>a.  (1) 

Here,  si,  is  the  maximum  refractive  index  of  the  core, 
A  is  the  relative  refractive  index  difference  between 
the  core  axis  and  cladding,  p  a  parameter  represent¬ 
ing  the  refractive  index  step  or  valley  at  the  core¬ 
cladding  boundary.  A  smooth  continuation  at  the 
core-dadding  boundary,  the  presence  of  a  step,  and 
that  of  a  valley  are  expressed  by  p- 1  ,p<  1,  andp>  1, 
respectively.  {a|aeR}  is  a  profile  parameter.  Some 
examples  of  a-index  profiles  are  shown  in  fig.  2. 

The  propagation  characteristics  of  an  optical  fiber 
are  governed  by  the  scalar  Helmholtz  differential 
equation  [IS] 


Fig.  I .  Optical  fiber  showing  grid  points  used  in  the  example. 
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Radius,  r 

Fit  2.  Examples  of  a-index  profiles,  a =oo  yields  a  step  index,  while  a = 2  yields  the  parabolic  index.  Values  of  p  control  the  character¬ 
istics  of  the  interface  between  the  core  and  the  cladding;  p<  1  results  in  a  step  at  the  interface,  p>  1  yields  a  valley. 


0+i£  + £)^0.  ,2, 

This  scalar  wave  equation  is  the  simplification  of  the 
exact  vector  wave  equation  under  the  assumption 
that  Vn/n  is  small,  which  includes  the  “small  index 
gradient’’  and  “weakly  guiding  approximations” 
[16,17].  In  the  above  equation  y/(r)  is  the  trans¬ 
verse  field  function  which  may  denote  either  the  di¬ 
electric  field  or  the  magnetic  field,  r  is  the  radial  co¬ 
ordinate,  n(r)  is  the  radial  refractive  index  profile, 
k  is  the  vacuum  wave  number,  fi  is  the  propagation 
constant  which  is  to  be  computed,  and  m  is  a  mode 
parameter  given  by 

m=  1  ,  for  TE and  TM  modes  (n=0) , 

=n+l,  for  EH  modes  (neN) , 

=  n—\,  for  HE  modes  (neN) .  (3) 

We  need  to  solve  the  differential  equation  in  order 
to  compute  the  propagation  constants.  From  the  ro¬ 
tational  properties  of  y  the  associated  boundary  con¬ 
dition  at  the  center  of  the  core  (r=0)  is 


(£L-°- form=0, 

y/(0)=0,  form^O.  (4) 

The  other  boundary  condition  applied  is  the  ex¬ 
tinction  of  field  at  the  jacket  written  as 

VW«  =  V'r_*=0,  (5) 

where  b  is  the  radius  of  core  and  cladding  together. 

2.1.  Transformation  to  nondimensional  form 


We  need  to  nondimensionalize  the  differential 
equation  for  easy  computation.  This  is  achieved  by 
setting 


u=y/y o,  x=r/a, 


(6) 


where  y0  is  the  maximum  field  amplitude  and  a  is 
the  radius  of  the  core.  Substituting  eq.  (6)  into  eq. 
(2)  we  obtain 
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By  including  the  refractive  index  distribution  given 
by  eq.  ( 1 ),  the  above  equation  can  be  rewritten  as 

SF  +  IS  +a‘  1  « 


where 

f(r)  =  2pd(r/a)a ,  OsSrsSa, 

=  2d  ,  a  <  r . 

Defining  the  parameters  U  and  W  as 
U=a(k2n2t  ~02) 1/2 ,  W=a(P2-k2n\)'n , 
we  can  define  V,  the  normalized  frequency,  as  [17] 
y2=U2+W2=k2a2(n2l-n22)  ,  (9) 

and  the  modified  propagation  constant,  #  as 

(10) 

Then  eq.  ( 8 )  becomes 

d2u  1  da  ^ j  V2n}  ,,  ,  m2\  „ 

s? + is + V~  ittAax)  -  "=0  ■ 


f(x)  =  2pAxa,  0«jc<1, 
=  2d  ,  jc>  1  . 


2.2.  Discretizing  the  differential  equation 

When  the  function  u  and  is  derivative  are  single 
valued,  finite  and  continuous  functions  of  x,  the  first 
and  the  second  differentials  can  be  approximated  by 
third  order  difference  formulas  as  follows  [18]: 

da  ui+ ,  —  , 

2*—-  (l3) 

d2«  u<+1  -2 _ 

—  ^  hi - •  04) 


«.  =  “(*) ,  Ui+x=u(x+h),  ut^l=u(x-h). 


h  is  the  width  between  the  grid  points  and  x=ih, 
{i  =  0,  1,  2, ...}.  Substituting eqs.  (13)  and  (14)  into 
(11),  and  defining 


«i -nj  ' 


we  get 

ui+x  -2Ui+Ui_i  1  Uf+i  —Ui_, 
h2  +  ih  2h 

+  0,  (16) 
and  on  rearranging,  the  equation  becomes 

■"-[-pO-s)] 

+W,[^(2+  7^~)  +  ^ lha )  “^] 
+““'[-p(1+3;)]“0'  ,l7) 

For  the  purpose  of  illustration,  we  have  chosen  six 
grid  points  along  the  radial  direction  as  shown  in  fig. 
1 .  In  general,  the  number  of  grid  points  can  be  any 
number  not  less  than  four,  the  minimum  necessary 
to  take  care  of  the  boundary  conditions.  Depending 
on  whether  m=0  or  m^O,  the  field  or  its  derivative 
vanishes  at  the  center  of  the  core.  When  the  deriv¬ 
ative  of  the  field  vanishes,  u0=ul. 

Writing  finite  difference  equations  at  the  grid 
points,  we  obtain  the  following  set  of  equations.  At 
i=l, 

“■  C+ ip — + *“>-*)  +u>{fi?)  -°  • 


where 

<5=  1  ,  m  =  0  , 

=  0,  m*  0. 


At  1  =  2, 


+  Pf(2ha) 
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At  i  =  3. 

=0. 

(20) 

At  i  =  4. 

(32  + m2  ,  r,rtAL  , 
k  l6h,  +  K/M *«>-J 

=0. 

(21) 

At  is  5,  since  the  field  goes  to  zero  at  the  jacket, 

(22) 

Finally,  at  is  6,  again  using  the  boundary  condition 
that  the  field  goes  to  zero  at  the  jacket, 

U6=  0.  (23) 

Since  the  boundary  condition  in  eq.  (23)  is  in¬ 
corporated  in  eq.  (22),  we  have  a  system  of  five 
equations. 


2.3.  Matrix  equation  formulation 


Formulating  a  matrix  equation  from  the  above  set 
of  equations  for  the  convenience  of  generalization 
and  easy  computation,  we  obtain 


/au-0  al2  \ 

fU'\ 

f  a2l  °21  —0  a23 

»2 

Aii  = 

<*32  a33  ~0  a3* 

Wj 

1  a**~0 

=0, 


(24) 


where  the  matrix  elements  are  defined  by 


a.j-, 

au  = 


—  (2/ —  1 ) 

» 

i  +  ?f(iha) , 


2/A 2 

2i2  +  2m2—5 


2  i2h 
2  i2  +  m2 
i2h 2 


+  Vf{  iha ) , 


i=  1  . 
i*  1  . 


ai.i+ 1  — 


-(2/-H) 
2  ih2 


(25) 


In  order  to  convert  the  problem  into  an  eigenvalue 
problem,  we  rewrite  eq.  ( 24 )  as 

( T-0l)u=Q ,  (26) 

where  I  is  the  identity  matrix,  and  T  is  a  tri-diagonal 
matrix. 

Equation  (26)  defines  an  eigenvalue  problem.  This 
means  that  eq.  (26)  has  a  nontrivial  solution  if  and 
only  if^is  an  eigenvalue  ( 19].  Hence,  the  required 
normalized  propagation  constants  contained  in  0 are 
obtained  by  finding  the  eigenvalues  of  the  tri-diag¬ 
onal  matrix  T.  This  mathematical  formulation  can 
be  generalized  to  {n  |  n  |  eN}  grid  points  in  the  radial 
direction  of  the  fiber  without  difficulty. 


3.  Arbitrary  profiles,  multiple  layers,  and  field 
distributions 

We  have  developed  this  method  of  analyzing  op¬ 
tical  fibers  using  the  a-index  profile.  This  is  because 
the  a-index  profile  is  commonly  used  in  the  litera¬ 
ture  and  can  represent  a  large  variety  of  real  refrac¬ 
tive  index  profiles,  including  the  very  important  step 
and  parabolic  profiles.  But  our  formulation  is  not 
limited  to  a-index  profiles. 

To  see  how  to  extend  the  method  to  arbitrary  pro¬ 
files  without  rederiving  a  system  of  finite  difference 
equations,  consider  eq.  (11).  The  refractive  index 
profile  is  included  in  this  equation  through  the  func¬ 
tion  f(x),  which  is  defined  in  eq.  ( 12).  Using/(x), 
the  refractive  index  profile  of  the  fiber  can  be  re¬ 
written  as 

n2(x)  =  n?(l-/(x)]  .  (27) 

Solving  for  f(x)  yields 

Rx)  =  \ -n2(x)/n2.  (28) 

By  generating  the  discretized  f=f(Xt=ih)  from 
samples  of  an  arbitrary  refractive  index  profile  n(x,), 
the  method  we  have  outlined  in  this  paper  can  be 
used  directly  on  arbitrary  profiles,  as  long  as  the 
“weakly  guiding”  approximation  holds. 

Multiple  layer  waveguides  of  any  number  of  layers 
may  be  considered  special  cases  of  arbitrary  refrac- 
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tive  index  profiles.  Since  we  have  normalized  the  fi¬ 
ber  core  radius,  a ,  to  unity,  we  must  explicitly  define 
which  layers  comprise  the  core  before  using  our 
method  and  scale  all  quantities  accordingly. 

Since  the  propagation  constant  for  a  mode  i,  /?„  is 
uniquely  associated  with  a  ft,  the  field  distributions 
for  a  propagating  mode  can  be  determined  from  the 
eigenvectors  «  of  matrix  T  (see  eq.  (26) ).  Many  ei¬ 
genvalue  routines  will  return  eigenvectors  as  well,  but 
at  the  cost  of  greatly  increasing  the  number  of 
computations. 

Useful  approximations  to  the  eigenvectors  for 
propagating  modes  can  be  computed  by  constructing 
the  tri-diagonal  matrix  T,  subtracting  a  specific  ft 
from  each  element  of  the  main  diagonal,  and  solving 
for  the  elements  of  it  using  standard  techniques  from 
linear  algebra.  From  the  observation  that  for  prop¬ 
agating  modes  the  field  will  approach  zero  at  the 
cladding/jacket  boundary,  we  can  set  uN,  the  right¬ 
most  element  of  u,  to  a  very  small  value  (not  zero), 
and  use  a  simple  backsubstitution  process  to  solve 
for  the  rest  of  the  u,.  This  procedure  yields  a  good 
approximation  to  the  field  distribution  multiplied  by 
an  arbitrary  constant. 


4.  Numerical  evaluation,  results,  and  discussion 

Although  the  mathematical  formulation  of  our 
method  for  determining  the  propagation  character¬ 
istics  of  an  optical  fiber  is  couched  in  terms  of  matrix 
equations,  there  are  special  structures  that  lead  to 
very  efficient  numerical  implementations.  First,  since 
T  is  a  tri-diagonal  matrix,  we  can  use  sparse  matrix 
techniques  to  reduce  storage  requirements  for  T. 
Second,  since  T  is  a  quasi-symmetric  tri-diagonal 
matrix,  we  can  use  a  similarity  transformation  to 
convert  T  into  a  real,  symmetric  matrix  [20].  Fi¬ 
nally,  the  eigenvalues  of  a  real,  symmetric  matrix  may 
be  computed  using  an  efficient  0(N2)  algorithm  (in 
our  case,  the  tqli.c  routine  from  ref.  [21  ],  which  has 
an  operation  count  of  approximately  30 N2). 

Using  eqs.  (23),  we  have  implemented  a  pair  of 
C  language  programs  which  compute  the  normalized 
propagation  constants  for  fibers  with  arbitrary  re¬ 
fractive  index  profiles.  We  define  the  normalized 
propagation  constant  as 
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U2  k.  n\  —  ft2 
X~  V 2  ~  k2n\—k2n\  ' 


(29) 


Note  that  some  authors  (e.g.  Gloge  [22] )  define  a 
normalized  propagation  constant  as 

b=\-U2/V2=\-x.  (30) 

One  program  computes  X  for  all  propagating  modes 
at  a  specific  value  of  m  in  eq.  (2)  over  a  range  of 
normalized  frequencies  V.  Another  program  searches 
for  the  cutoff  frequency  (Kc)  of  a  specific  linearly 
polarized  ( LP )  mode.  Both  programs  allow  all  of  the 
parameters  in  eq.  ( 1 )  to  be  varied,  as  well  as  the  val¬ 
ues  of  b  and  N„  which  are  the  fiber  radius  (see  fig. 
1)  and  number  of  grid  points  in  the  fiber  core, 
respectively. 

In  verifying  the  performance  of  our  method,  we 
have  computed  the  propagation  characteristics  of  step 
index  and  parabolic  index  fibers  over  a  normalized 
frequency  range  of  0  to  20.  These  index  profiles  have 
analytical  solutions  and  have  been  studied  analyti¬ 
cally  and  numerically  by  other  authors  [  14,23-26]. 
Our  results  agree  well  with  previously  published  re¬ 
sults,  as  shown  in  table  1.  Note  that  for  propagating 
modes,  x  must  lie  between  0  and  1  (i.e.  0^*$  1 ). 

For  comparison  with  a  known  case,  fig.  3  shows 
the  dispersion  characteristics  we  have  computed  for 
the  step  index  profile.  The  plot  agrees  well  with  the 
analytic  results  for  fibers  with  infinite  cladding.  The 
small  differences  between  the  computed  and  analytic 
cutoff  frequencies  for  each  mode  are  due  to  the  finite 
cladding  width  used  in  our  computations,  and  the 
finite  number  of  grid  points  across  the  fiber.  The 
fundamental  mode,  which  has  zero  cutoff  in  the  in¬ 
finite  cladding  case,  is  extremely  sensitive  to  the 
cladding  width. 

Figures  4,  5,  and  6  show  our  dispersion  calcula¬ 
tions  for  the  parabolic  index  fiber  with  cladding  width 
10  times  the  core  width  and  p=  1.0,  2.0,  and  0.75. 
These  compare  favorably  with  published  results  (see, 
for  example,  refs.  [  14,23,26] ). 

Two  factors  have  a  major  influence  on  the  results 
of  our  computations:  the  number  of  grid  points  used 
across  the  fiber  (which  we  specify  in  terms  of  the 
number  of  points  in  the  core  of  the  fiber),  and  the 
width  of  the  cladding.  Figure  7  shows  the  results  of 
applying  our  method  to  a  step  index  fiber  for  two  LP 
modes.  For  each  mode  we  have  calculated  the  cutoff 
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Table  1 

Comparison  of  the  cutoff  frequencies  obtained  by  the  finite  difference  method  with  analytical  and  previous  numerical  results.  R is 
the  fiber  radius,  R ^  is  the  core  radius,  and  <5  is  the  percentage  difference  from  the  infinite  cladding  result.  256  points  were  used  in  the 
fiber  core. 


a 

Mode  (m,  1) 

Infinite  cladding 

Normalized  cutoff  frequency 

^  fiber  =  10-Rcore 

<5(%) 

^  fiber  — 

<S(%) 

i 

1,1 

4.381 

4.391 

0.23 

4.384 

0.07 

2 

1,1 

3.518 

3.526 

0.23 

3.520 

0.06 

1,2 

7.451 

7.457 

0.08 

7.453 

0.03 

2,1 

5.744 

5.744 

<10'2 

5.744 

<10"2 

2,2 

9.645 

9.645 

<10“2 

9.645 

<  10~2 

3,1 

7.848 

7.848 

<10“2 

7.848 

<  10~2 

4,1 

9.904 

9.904 

<10'2 

9.904 

<10-2 

3 

1,1 

3.181 

3.189 

0.3 

3.183 

0.06 

4 

1.1 

3.000 

3.007 

0.2 

3.002 

0.07 

5 

1.1 

2.886 

2.894 

0.28 

2.888 

0.07 

10 

1.1 

2.649 

2.657 

0.30 

2.651 

0.08 

20 

1,1 

2.527 

2.535 

0.32 

2.529 

0.08 

oo 

1.1 

2.405 

2.413 

0.33 

2.407 

0.08 

Fig.  3.  Dispersion  characteristics  of  a  step  index  fiber  (a=oo,  0.038). 


frequency  using  from  4  to  256  points  in  the  core,  and 
for  fiber  radii  from  5  to  20  times  the  core  radius. 
From  this  figure  we  can  see  the  expected  conver¬ 
gence  on  a  final  result  as  the  number  of  points  in  the 


core  increases.  The  effect  of  cladding  width  is  also 
apparent. 

The  effect  of  the  number  of  grid  points  in  the  core 
is  two  fold.  As  the  number  of  grid  points  is  increased 
the  distance  between  samples  of  the  refractive  index 
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profile  is  decreased,  resulting  in  a  better  approxi¬ 
mation  of  the  actual  profile.  This  is  especially  ap¬ 
parent  in  profiles  with  sharp  transitions  at  the  core¬ 
cladding  interface,  such  as  for  a »  I.  Cases  where 
p^\  (see  eq.  ( 1 ) )  are  also  likely  to  be  poorly  mod¬ 
eled  by  a  small  number  of  core  sample  points.  The 
effect  of  reducing  the  number  of  sample  points  in  the 
core  is  to  apply  a  “low  pass”  spatial  filter  to  the  re¬ 
fractive  index  profile. 

Setting  the  number  of  points  in  the  core  also  ef¬ 
fectively  applies  a  filter  to  the  spatial  frequency  con¬ 
tent  of  the  field  distributions  calculated  for  each 
mode  in  the  fiber.  When  computing  propagation 
constants  at  higher  normalized  frequencies,  using  a 
small  number  of  samples  may  induce  errors  due  to 
a  form  of  “aliasing”.  These  two  effects  are  respon¬ 
sible  for  the  poor  results  when  the  number  of  grid 
points  in  the  core  is  below  approximately  10  for  the 
modes  we  have  examined. 

Using  the  step  index  fiber  as  an  example,  fig.  8 
demonstrates  the  effects  of  the  number  of  grid  points 
by  plotting  the  computed  cutoff  frequency  for  modes 
LP01,  LP,,,  and  LP02  for  several  different  grid  sizes. 
In  this  figure,  each  curve  is  normalized  to  the  value 


of  the  cutoff  frequency  for  that  mode  calculated  with 
256  points  in  the  core.  We  can  see  that  for  mode  LP01 
the  cutoff  frequency  calculated  with  8  points  in  the 
core  is  less  than  1.0015  times  that  computed  using 
256  points  in  the  core,  while  for  mode  LP,,  (with  a 
higher  cutoff  frequency)  we  need  at  least  12  points 
in  the  core  for  similar  results.  In  general,  as  the  nor¬ 
malized  frequency  increases,  the  number  of  points  in 
the  core  must  be  increased  to  maintain  the  accuracy 
of  the  method. 

For  modes  with  relatively  low  cutoff  frequencies, 
variations  in  cladding  width  produce  large  changes 
in  the  calculated  cutoff  frequency,  Cutoff  fre¬ 
quency  increases  as  the  cladding  width  decreases. 
This  is  the  expected  behavior.  Analyses  assuming  in¬ 
finite  cladding  width,  while  adequate  for  many  pur¬ 
poses,  fail  to  account  for  the  increasing  importance 
of  finite  cladding  width  as  the  normalized  frequency 
becomes  smaller.  The  fundamental  mode,  which  has 
no  cutoff  frequency  when  the  cladding  is  infinite, 
shows  a  definite  cutoff  in  real  fiber. 

Figure  9  shows  the  effects  of  cladding  width  on  the 
cutoff  frequencies  of  two  LP  modes  in  a  step  index 
fiber.  In  this  plot,  the  curves  for  each  mode  are  nor- 


Fig.  8.  Effect  of  the  number  of  grid  points  on  the  computed  cutoff  frequencies  of  propagating  modes  of  a  step  index  fiber  (LP0i.  LP,,, 
and  LP02). 


402 


132 


Volume  99,  number  5,6 


OPTICS  COMMUNICATIONS 


15  June  1993 


malized  to  the  cutoff  frequency  for  that  mode  in  the 
infinite  cladding  case.  As  an  example,  consider  mode 
LP,,.  We  can  see  that  the  cutoff  frequency  when  the 
fiber  radius  is  10  times  the  core  radius  is  only  about 
1.003  times  (0.3%)  greater  than  the  cutoff  frequency 
when  the  fiber  radius  is  25  times  the  core  radius. 
However,  when  the  fiber  radius  is  only  5  times  the 
core  radius  then  the  cutoff  frequency  increases  to 
1 .0 1 4  times  (1.4%)  greater  than  the  cutoff  frequency 
when  the  fiber  radius  is  25  times  the  core  radius. 

Using  the  approach  outlined  in  sect.  3,  we  have 
computed  the  field  distributions  for  three  LP  modes 
in  a  step  index  fiber.  Figure  10  shows  the  results,  with 
the  field  patterns  normalized  so  that  the  maximum 
value  in  each  pattern  is  one.  These  results  agree  well 
with  the  results  reported  in  ref.  (15). 

5.  Conclusions 

We  know  that  when  A  is  small  in  an  optical  fiber, 
the  scalar  approximation  yields  results  that  are  very 
close  to  the  exact  vector  formulation.  Even  for  large 
differences  between  the  core  and  cladding  refractive 
indexes,  optimum  single-mode  fiber  parameters  ob¬ 
tained  from  the  scalar  approximation  differ  negli¬ 
gibly  from  those  obtained  using  the  exact  formula¬ 
tion  [27]. 

We  have  developed  a  method  to  evaluate  the  prop¬ 
agation  constants  by  transforming  the  scalar  wave 
equation  into  a  set  of  finite  difference  equations  and 
then  converting  into  a  matrix  eigenvalue  problem. 
The  method  does  not  involve  a  search  procedure  to 
find  the  propagation  constants,  or  the  explicit  eval¬ 
uation  of  Bessel  and  modified  Bessel  functions,  which 
is  time  consuming,  as  was  the  case  in  earlier  works. 

We  have  demonstrated  the  convergence  of  the 
method  and  the  dependence  of  the  rate  of  conver¬ 
gence  on  the  number  of  grid  points.  The  method  is 
elegant,  stable,  straight  forward,  is  applicable  to  ar¬ 
bitrary  index  profiles,  and  is  accurate.  We  have  also 
explored  and  established  the  effects  of  finite  clad¬ 
ding  width  on  the  dispersion  characteristics  of  op¬ 
tical  fiber. 
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ABSTRACT 

We  have  developed  a  set  of  computer  codes  that  compute  the  propagation  constants 
and  field  patterns  for  the  propagating  modes  of  cylindrical  optical  fibers.  From  a  simple 
set  of  finite  difference  equations,  solutions  of  the  scalar  Helmholtz  wave-equation  may 
be  computed  across  a  range  of  normalized  frequencies  to  generate  curves  describing 
the  dispersion  characteristics  of  the  fiber.  Accurate  cutoff  frequencies  for  any  mode  can 
also  be  computed.  We  designed  the  computer  codes  around  a-index  profiles  since  these 
profiles  have  been  extensively  covered  in  the  literature,  but  our  system  also  supports 
arbitrary  profiles  within  the  limits  of  the  "small  index  gradient”  and  "weakly  guiding" 
approximations.  The  computer  codes  are  accurate  and  fast.  They  may  be  used  interactively 
to  explore  dispersion  in  optical  fibers  and  the  effects  of  finite  cladding  width  on  dispersion. 
©  1993  |ohn  Wiley  &  Sons.  Inc. 


INTRODUCTION 

By  simplifying  a  set  of  computer  codes  that  we  use 
in  our  research,  we  have  put  together  a  system  of 
computer  codes  for  use  in  illustrating  the  propaga¬ 
tion  characteristics  of  optical  fibers  to  both  graduate 
students  and  undergraduates.  Comprising  a  menu 
interface,  computational  kernel,  and  graphical  dis¬ 
play,  our  pedagogical  system  can  run  on  most  IBM- 
PC  compatible  computers  with  a  graphics  adaptor 
and  numerical  coprocessor. 


MATHEMATICAL  FORMULATION 

Optical  fibers  are  inhomogeneous  dielectric  cylin¬ 
ders  of  radius  a,  called  the  “core,”  surrounded  by 
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a  homogeneous  refractive  index  medium  called  the 
“cladding.”  The  cladding,  in  turn,  is  encased  in  a 
highly  lossy  material  called  the  “jacket.” 

The  refractive  index  profile  of  the  fiber,  called  an 
a-index  profile,  is  given  by 

\n7,[\  -  2pA(r/a)a]  for  Osrsa 

n V)  =  {  , 

[n?[l  -  2A]  for  r>  a 

(1) 

Here,  nt  is  the  maximum  refractive  index  of  the 
core,  A  is  the  relative  refractive  index  difference  be¬ 
tween  the  core  axis  and  cladding,  and  p  is  a  param¬ 
eter  representing  the  refractive  index  step  or  valley 
at  the  core-cladding  boundary.  A  smooth  contin¬ 
uation  at  the  core-cladding  boundary  is  expressed 
by  p  =  I,  while  the  presence  of  a  step  is  indicated 
by  p  <  I .  Setting  p  >  I  results  in  a  valley  at  the  core¬ 
cladding  interface.  Some  examples  of  a-index  pro¬ 
files  are  in  Figure  I . 
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Figure  I  Examples  of  a-index  profiles. 


The  propagation  characteristics  of  an  optical  fiber 
are  governed  by  the  scalar  Helmholtz  differential 
equation  [1] 


dr 2  +  ~r  dr  + 


n2(r)k 2  -  P2  - 


(2) 


This  scalar  wave  equation  is  the  simplification  of 
the  exact  vector  wave  equation  under  the  “small 
index  gradient”  and  “weakly  guiding”  approxima¬ 
tions  [1J.  When  A  1  in  an  optical  fiber,  the  scalar 
approximation  yields  results  that  are  very  close  to 
the  exact  vector  formulation,  and  even  for  larger 
differences  between  the  core  and  cladding  refractive 
indexes,  optimum  single-mode  fiber  parameters  ob¬ 
tained  from  the  scalar  approximation  differ  negli¬ 
gibly  from  those  obtained  using  the  exact  formu¬ 
lation  [2],  For  fibers  used  in  communication  ap¬ 
plications,  A  <  0.03  is  common  [3], 

In  Eq.  (2),  \fir)  is  the  transverse  field  function 
(eiiher  the  electric  field  or  the  magnetic  field),  r  is 
the  radial  coordinate,  n(r)  is  the  radial  refractive 
index  profile,  k  is  the  vacuum  wave  number,  /3  is 
the  propagation  constant  which  is  to  be  computed, 
and  m  is  a  mode  parameter  given  by 


f  1  for  TEo/  and  TM0/  modes  (p  =  0) 


m  = 


p  +  1  for  EHp/  modes  (p  G  N) 


-  1  for  HEP,  modes  (p  G  N) 


(3) 


Here,  N  is  the  set  of  natural  numbers. 

Modes  with  the  same  propagation  constants  are 
grouped  under  a  linearly  polarized  ( LP)  mode  clas¬ 
sification.  Each  propagating  mode  is  identified  by 
an  LPm,  designation  where  m  is  defined  in  Eq.  (3) 
and  /  is  identical  to  the  value  in  the  traditional  HEW, 
EHpj,  TEp/,  and  TM^  mode  designation  [3). 

We  need  to  solve  the  differential  equation  to 
compute  the  propagation  constants.  From  the  ro¬ 
tational  properties  of  \p  the  associated  boundary 
condition  at  the  center  of  the  core  (r  =  0)  is 


=  0  for  m  =  Q 

(4) 

^(0)  =  0  for  m  +  0 

The  other  boundary  condition  applied  is  the  ex¬ 
tinction  of  field  at  the  jacket,  written  as 

^jack«  =  ypr.b  =  o  (5) 

where  b  is  the  radius  of  core  and  cladding  together. 

The  issue  of  boundary  conditions  is  complex  but 
very  important  in  all  numerical  work.  A  more  ap¬ 
propriate  boundary  condition  in  the  jacket  is  an  ab¬ 
sorbing  boundary  condition,  but  for  a  large  cladding 
width,  as  assumed  here,  our  boundary  condition  is 
appropriate  [I],  For  an  excellent  discussion  of  ra¬ 
diation  boundary  conditions,  see  Moore  et  al.  [4], 


dip 
dr  . 
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We  need  to  rewrite  the  differential  equation  to 
eliminate  dependence  on  units;  this  is  a  more  general 
formulation  that  is  also  easier  for  computation.  We 
do  this  by  setting 


with 


/(*)  = 


2pA.\" 

2A 


0  <  A 1 
X  >  1 


(12) 


and 


When  the  function  u  and  its  derivative  are  single 
valued,  finite  and  continuous  functions  of  x,  the 
first  and  the  second  differentials  can  be  approxi¬ 
mated  by  second  order  difference  formulas.  Using 
the  finite  difference  approximation,  and  defining 


x  = 


r 

a 


(6) 


V  = 


V2n] 

(rti  ~  n\) 


(13) 


where  4>o  is  the  maximum  field  amplitude  and  a  is 
the  radius  of  the  core.  Substituting  Eq.  (6)  into  Eq. 
(2)  we  obtain 


d2u  1  du 
dx2  +  x  dx 

+  a2j/c2/i2(xa)  -  02  -  =  0  (7) 

By  including  the  refractive  index  distribution 
given  by  Eq.  ( 1 ),  and  normalizing  all  waveguide 
dimensions  to  the  core  radius,  the  above  equation 
can  be  rewritten  as 

d2u  |  du~ 
dx2  +  x  dx 

+  [*2/i?(l  -/(*)) -/32-^ju  =  0  (8) 

Defining,  V,  the  normalized  frequency,  as 

V2  =  k2[n]  -  n22]  (9) 

and  the  modified  propagation  constant,  0,  as 


0  = 


V2n2 


^  “ 


@2  (10) 


we  arrive  at  the  discretized  wave  equation 


(14) 


where  h  is  the  distance  between  grid  points  and  x 
=  ih,  {/  =  0,  1,  2*  •  • }. 

Writing  finite  difference  equations  at  the  grid 
points,  we  obtain  a  set  of  equations  that  may  be 
written  as  a  matrix  equation: 


Au  =  0  (15) 

To  convert  the  problem  into  an  eigenvalue  problem, 
we  rewrite  Eq.  ( 15)  as 


[T-/?I]u  =  0  (16) 

where  I  is  the  identity  matrix,  and  T  is  a  tri-diagonal 
matrix.  Equation  (16)  has  a  nontrivial  solution  if 
and  only  if  0  is  an  eigenvalue  [  5  ] .  Hence,  the  re¬ 
quired  normalized  propagation  constants  contained 
in  0  are  obtained  by  finding  the  eigenvalues  of  the 
tridiagonal  matrix  T. 


Eq.  ( 8 )  becomes 


EXAMPLES 


d2u  1  du 
dx2  +  x  dx 


0~ 


V2ni 


(ni  ~  n]) 


f(x) 


m 


■  u  =  0 


(ID 


To  construct  a  system  of  computer  codes  that  will 
run  well  on  IBM-PC  compatible  computers,  we  have 
taken  advantage  of  the  special  properties  of  our  ma¬ 
trix  formulation.  Most  importantly,  since  T  is  a 
quasi-sym metric  tri-diagonal  matrix,  we  can  use  a 
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Version  2.14 

Fiber  Parameter  Selection 

Numerical  Procedures 

Profile  Type: 

Computation: 

Step: 

X 

Dispersion: 

X 

Graded: 

Field  Pattern: 

Arbitrary: 

Cutoff  Frequency: 

Profile  parameters: 

Computational  Parameters: 

Alpha : 

Starting  V: 

00-2S 

Delta: 

tuna 

Ending  Vs 

20.00 

Rho : 

Step  Site  for  V: 

00.25 

Core  Refractive  Idx: 

1-S20 

Mode  Number  (m,  1) : 

00 

Clad  Refractive  Idx: 

Cutoff  Lower  Bound: 

Core  Radius  (um) : 

01. 0Q 

Cutoff  Upper  Bound: 

Fiber  Radius  (um) : 

i£Ljm 

Points  in  Core: 

022 

Figure  2  Examples  of  the  Fiber  User  Interface  for  a  step  index  fiber. 


similarity  transformation  to  convert  T  into  a  real, 
symmetric  matrix  [6].  The  eigenvalues  of  a  real, 
symmetric  matrix  may  be  computed  using  an  effi¬ 
cient  0(N2)  algorithm  (in  our  case,  the  tql  i  rou¬ 
tine  from  [7],  which  has  an  operation  count  of  ap¬ 
proximately  30  N2). 


The  computational  kernel  of  our  pedagogical 
system  consists  of  a  trio  of  computer  codes  that 
compute  the  normalized  propagation  constants  and 
field  patterns  for  cylindrical  fibers.  Arbitrary  refrac¬ 
tive  index  profiles  (which  must  meet  the  criteria  for 
Eq.  (2)  to  be  valid)  are  read  from  a  file  on  disk, 


Figure  3  Normalized  propagation  constant  vs.  normalized  frequency  for  a  step  index  fiber  ( a 
=  oo,  A  =  0.038).  =  10  R<on  and  =  32. 
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Figure  4  Field  patterns  for  the  modes  LP0i ,  LPu,  and  LP2,  of  a  step  index  fiber  (or  =  oo,  A 
=  0.038)  at  V  =  10.  /W  =  10K««  and  Nc  =  32. 


while  a-index  profiles  are  constructed  directly  within 
each  program.  The  numerical  computer  codes  call 
upon  a  common  graphical  display  program  and  file 
I/O  routines.  A  single  shell  program  controls  the 
entire  suite  of  computer  codes  and  coordinates  the 
loading  and  execution  of  code  segments  as  required 
to  perform  the  computations  requested  by  a  user. 

Common  to  all  computer  codes  in  the  suite,  the 
normalized  propagation  constant  is  defined  by 


For  propagating  modes,  kn2  Zen,  [3],  and  so 

X  must  lie  between  0  and  I.  All  computer  codes 
allow  the  parameters  in  Eq.  ( I )  to  be  varied,  as  well 
as  the  values  of  b  and  Nc,  the  fiber  radius  and  num¬ 
ber  of  grid  points  in  the  fiber  core,  respectively. 

To  illustrate  our  system,  we  have  computed  the 
propagation  characteristics  of  a  step  index  (a  =  oo) 
and  a  parabolic  index  fiber  (a  =  2)  over  a  normal¬ 
ized  frequency  range  of  0-20.  These  refractive  index 
profiles  have  solutions  that  have  been  studied  ana¬ 
lytically  and  numerically  by  other  authors 
18,9,10,1 1,12),  and  our  results  agree  well  with  pre¬ 
viously  published  results. 

The  first  step  in  an  analysis  is  to  define  the  fiber 
profile  parameters,  and  the  number  of  points  to  use 
in  finite  difference  approximations.  These  defini¬ 
tions  are  made  in  the  shell,  as  shown  in  Figure  2  for 


the  step  index  case.  Then,  a  user  indicates  which 
computation  is  desired.  The  shell  runs  appropriate 
computer  codes  to  generate  the  desired  data  and 
graphs. 

Figure  3  shows  the  step  index  fiber  dispersion 
curves  generated  using  Nc  -  32,  and  Figure  4  shows 
the  field  patterns  for  the  LP0, ,  LP, , ,  and  LP2,  modes. 
The  fiber  parameters  are  as  shown  in  Figure  2.  Field 
intensity  patterns  for  each  LP  mode  are  defined  by 


u 2 

where  um„  is  the  maximum  field  magnitude  along 
the  fiber  radius  for  each  LP  mode. 

Dispersion  curves  for  a  parabolic  index  fiber  are 
shown  in  Figure  5.  These  curves  were  generated  us¬ 
ing  a  =  2  and  p  =  1;  all  other  parameters  were  iden¬ 
tical  to  those  for  the  step  index  case  presented  in 
Figure  3.  Figure  6  shows  the  field  pattern  for  the 
LP0| ,  LP,  i ,  and  LP2|  modes. 

Figures  3-6  were  prepared  using  data  passed  di¬ 
rectly  to  gnuplot,  a  powerful  scientific  function 
and  data  graphing  program  available  without  charge 
on  a  variety  of  platforms,  including  IBM-PC  com¬ 
patibles.  By  using  gnuplot  as  our  graphical  display 
routine,  we  allow  students  to  view  plots  on  any  IBM- 
PC  with  a  standard  graphics  adapter  and  prepare 
publication  quality  graphics  on  any  of  the  devices 
supported  by  gnuplot. 
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Figure  S  Normalized  propagation  constant  vs.  normalized  frequency  for  a  parabolic  index 
fiber  (a  =  2.  p  =  1,  A  =  0.038).  Rtbet  =  10/?^  and  Nc  =  32. 


Using  our  suite  of  computer  codes,  students  can 
compute  dispersion  characteristics  and  field  patterns 
for  a  variety  of  refractive  index  profiles  in  one  ses¬ 
sion.  The  computations  and  graphs  can  be  printed 
and  compared,  facilitating  an  understanding  of  the 
effect  on  propagation  of  varying  fiber  parameters. 


Due  to  the  nature  of  solutions  that  can  be  com¬ 
puted  directly  from  the  scalar  Helmholtz  wave 
equation,  the  field  patterns  generated  by  our  system 
are  limited  to  showing  the  radial  variations  of  each 
mode.  Students  may  determine  the  angular  variation 
of  a  field  pattern  by  relating  the  LPm/  mode  desig- 


Figure  6  Field  patterns  for  modes  LPoi ,  LP,, .  and  LPji  of  a  parabolic  index  fiber  (a  =  oo ,  p 
=  I,  A  =  0  038)  at  V  =  10.  R ^  =  1 0Rcon  and  N<  =  32. 
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nations  to  TE,  TM,  HE,  and  EH  mode  designations 
as  described  in  the  introduction,  and  noting  that 
there  are  2m  field  maxima  around  the  fiber  circum¬ 
ference  and  /  field  maxima  along  the  fiber  radius 
[  3  ] .  A  further  release  of  our  software  will  include 
modules  designed  to  aid  in  the  construction  of  full 
2-D  field  patterns. 

Two  factors  have  a  major  influence  on  the  results 
of  our  computations:  the  number  of  grid  points  used 
across  the  fiber  (which  we  specify  as  the  number  of 
points  in  the  core  of  the  fiber),  and  the  width  of  the 
cladding.  For  modes  with  low  cutoff"  frequencies, 
variations  in  cladding  width  produce  large  changes 
in  the  calculated  cutoff"  frequency,  Vc.  Cutoff"  fre¬ 
quency  increases  as  the  cladding  width  decreases. 
This  is  the  expected  behavior.  The  fundamental 
mode,  which  has  no  cutoff  frequency  when  the 
cladding  is  infinite,  shows  a  definite  cutoff  in  real 
fiber. 

The  number  of  grid  points  across  the  fiber  affects 
the  accuracy  of  the  finite  difference  approximations 
used  in  computing  solutions  to  the  wave  equation. 
In  general,  a  coarse  grid  results  in  an  apparent  shift 
of  all  propagation  constants  that  increases  with  fre¬ 
quency.  In  our  experience,  16  points  in  the  core  are 
sufficient  for  normalized  frequencies  below  10.  Nc 
-  1 28  is  sufficient  for  most  cases  of  interest,  but  will 
run  slowly  even  on  a  ’486  based  computer. 

CONCLUSIONS 

We  have  developed  an  integrated  set  of  computer 
codes  to  evaluate  propagation  constants  and  field 
patterns  of  modes  by  transforming  the  scalar  wave 
equation  into  a  set  of  finite  difference  equations  and 
then  converting  into  a  matrix  eigenvalue  problem. 
Our  computer  codes  are  fast  enough  to  run  on  an 
IBM-PC  with  a  numeric  coprocessor,  are  accurate, 
and  provide  a  convenient  system  with  which  stu¬ 
dents  can  explore  propagation  in  optical  fibers. 
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